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ABSTRACT 

In this work, a low-noise plasma simulation model is developed, and 
applied to a series of linear and nonlinear problems associated with 
electrostatic wave propagation in a one-dimensional, collisionless, 
Maxwellian plasma, in the absence of magnetic field. It is demonstrated 
that use of the hybrid simulation model allows economical studies to be 
carried out in both the linear and nonlinear regimes with better quanti- 
tative results, for comparable computing time, than can be obtained by 
conventional particle simulation models, or direct solution of the 
Vlasov equation. 

The characteristics of the hybrid simulation model itself are first 
investigated, and it is shown to be capable of verifying the theoretical 
linear dispersion relation at wave energy levels as low as 10~® of the 
plasma thermal energy. Having established the validity of the hybrid 
simulation model, it is then used to study the nonlinear dynamics of a 
monochromatic wave, sideband instability due to trapped particles, and 
satellite growth. The simulations are performed in parameter ranges 
such that detailed quantitative comparison with available theories is 
possible. In particular, the transition from time-asymptotic amplitude 
oscillation to continuous Landau damping is investigated for a mono- 
chromatic wave as the initial wave amplitude is varied. The results, 
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which include a small nonlinear frequency shift, compare favorably with 
theory. The study of sideband instability confirms the applicability 
of quasilinear theory at short times, and parametric coupling theory in 
the time-asymptotic limit, and reproduces features analogous to those 
observed in laboratory experiments. The growth of a satellite wave in 
the presence of a large amplitude wave and a test wave is shown to be 
explicable by a simple theory involving slow modulation of the large 
amplitude wave. 

It is concluded that the hybrid simulation model constitutes a 
reliable and economical tool for use in the study of plasma phenomena. 

It should be widely applicable to the testing of predictions of linear 
and nonlinear theories, and the description of more complicated situations 
which are not amenable to analysis. 
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1. INTRODUCTION 


This dissertation is concerned with computer simulation of electron 
plasma waves in a Maxwellian plasma, in the absence of magnetic field. 
The aims of our simulations are, first, to develop an economical low- 
noise simulation technique, and second, to apply it to the study of 
linear and nonlinear wave phenomena in a one-dimensional plasma with 
periodic boundary conditions. In what follows, emphasis has been placed 
on simulations for wave and plasma parameters comparable with those 
assumed in available theories, and accessible in laboratory experiments. 

There are two distinctly different approaches to the simulation of 
plasma dynamics: first is the use of a particle simulation model in 
which individual charged particles are followed, and second is direct 
numerical solution of the Vlasov equation describing the charged particle 
velocity distribution function. The particle simulation model has the 
disadvantage that the fluctuation level is usually several orders of 
magnitude higher than in an actual plasma. This stems from the fact 
that it is not feasible to follow on the computer the dynamics of as 
many particles as there are in a plasma. This fluctuation not only gives 
rise to nonphysical effects, but also makes it difficult to study linear 
and weakly nonlinear phenomena. This is particularly unfortunate since 
most of the nonlinear theories to date are based on an expansion method 
which is valid only in weakly nonlinear cases. Consequently, they 
cannot be clearly validated by computer simulation, nor vice versa. 

Direct solution of the Vlasov equation is subject to numerical instabi- 
lity associated with the free-streaming term in the Vlasov equation. 

This tends to limit application of the method to short-time simulations. 
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It will be shown in this work that a hybrid approach, combining 
features of the particle simulation model and Vlasov approach, can 
avoid the foregoing difficulties, and provide reliable results in both 
the linear and nonlinear regimes. 

In Section 2, we review various types of particle simulation models 
and consider numerical solution of the Vlasov equation. The hybrid 
simulation model is described and analyzed in Section 3, and results of 
test runs in the linear regime are presented to demonstrate the feasi- 
bility of high quality simulation. In Section 4, the nonlinear dynamics 
of a monochromatic plasma wave are studied as our first application of 
the model. Section 5 presents a comprehensive study of sideband insta- 
bility in which a large amplitude wave and its sidebands interact 
through the trapped particles to cause growth of the sideband waves. 

In Section 6, satellite growth due to nonlinear interaction between the 
lower and upper sidebands is studied. Some conclusions on the validity 
and applicability of the hybrid simulation approach are given in 
Section 7* 
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2. COMPUTER APPLICATIONS TO PLASMA PHYSICS 
2 • 1 The Vlasov Equation and Simulation Approaches 

There are two ways of approaching problems in plasma physics by 
use of computers. One is to use the Vlasov equation, describing the 
velocity distribution function of smeared-out charges, combined with 
the Maxwell equations. The other is to use a particle simulation 

approach dealing with discrete charges, combined with the Newton equation 
and the Maxwell equations . 

The Vlasov equation, 


a± 'i . a*. 

& r + z ‘fc" + ^~ + ~ x £ > ’ -^r = 0 > (2.D 


coupled with the Maxwell equations, 


^ x £ = £ + TT > V • B = 0 , 

at 

( 2 . 2 ) 

SB 

7 x £- -st ■ 7 • E-f- • 

0 


and equations for current and charge densities, 


1 ~2/' q i~ f i d ~ ’ P = E/*i f i d ~ * <2. 3) 

i i 

is capable of describing the collective behavior of collisionless 

plasmas. Here, q is the distribution function of the i-th charged 

particle species, E and B are the electric and magnetic fields, 

including both time-varying and static fields, q. and m are the 

i i 

particle charge and mass, and e Q is the permittivity of free space. 
It is important to note that the Vlasov equation does not include 
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discrete particle encounters and thermal fluctuations, whereas discrete 

charge models automatically include such effects. 

In the Vlasov approach, the major task is to solve Eqs. (2,l)-(2.3) 

numerically, subject to appropriate initial and boundary conditions. 
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This line has been pursued by Knorr, Kellogg, and Armstrong. These 
authors solved the Vlasov equation, 


&f. &f. q. af. 

i i 

+ v - — + — E -r — 

at ax in, av 


0 , 


(2.4) 


and the Poisson equation, 


dE 

dx 



eA'i*' • 


(2.5) 


1 3 

in the study of Landau damping, y and the two-stream instability at 
2 

large amplitudes, in a one-dimensional plasma. However, all of these 
workers encountered the serious computational difficulty which stems 
from the nature of the Vlasov equation, i.e,, the tendency to develop 
steep gradients in phase-space as time increases. To appreciate this, 
consider the following simplified equation, 


3f 

at + 



( 2 . 6 ) 


which describes free- streaming particles. By taking a single Fourier 
component of f , it is seen from Eq. (2.6) that the distribution 
function always has a part of the form, exp (ikvt). This represents 
velocity- space oscillation which becomes finer and finer as time 
increases. It is due to conversion of the space oscillations of the 
initial distribution function into velocity oscillations by the shear 
motion in phase-space. 
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Since the gradient of the velocity distribution function increases 

with time, a progressively finer grid in phase-space will be required 

for accurate computation. This is the case for the difference method 

used by Kellogg. in the case of the expansion method used by Knorr 1 
3 

and Armstrong, the development of large derivatives of f requires a 
large number of expansion terms to be used. The number of grid points, 
or the number of expansion terms allowable, is limited by computer 
capacity and economic considerations. The finite grid size, or trun- 
cated expansion, consequently sets limits on the elapsed time during 
which the simulation is accurate. 

Direct solution of the Vlasov equation gives accurate quantitative 
results in both the linear and nonlinear regimes only until it approaches 
the limit mentioned above. The particle simulation model, on the other 
hand, is free from such time limitations, and is therefore better 
suited for nonlinear problems. 

Computer simulation seems to have been introduced into the field 
of plasma physics by Buneman, who used the charge sheet model to solve a 
nonlinear problem . 4 In particular, he followed the development of 
electron-ion two-stream instability to show that rapid randomization of 
the initially coherent streaming energy takes place in a cold collision- 
less plasma. Soon afterwards, Dawson investigated the thermalization 
and ergodic behavior of a system of charge sheets and found that 
statistical mechanics can be applied to the system . 5 He then studied 
the properties of the one-dimensional plasma in thermal equilibrium, 
finding good agreement with theoretical results based on statistical 
mechanics. Since then, the simulation model has been improved and 
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applied to a wide variety of plasma problems which are not analytically 
tractable. Computer codes have been written for both electrostatic 
and electromagnetic cases, with or without external magnetic field, 
and in one, two, or three dimensions* Relativistic effects have also 
been included. 

In the charge sheet model, electrons or ions are considered to 
form infinitesimally thin charge sheets, extending to infinity in two 
dimensions, and moving only in the direction perpendicular to the plane 
of the sheet. They are assumed to pass freely through each other. Such 
a charge sheet generates electric field, E , which is discontinuous at 
the sheet position, and takes constant values on each side, as shown in 
Fig. 2.1. For a given boundary value, the electric field created by 
many charge sheets can be calculated at any point by a simple algorithm. 
For example, suppose the electric field to be given at x = 0 . The 
algorithm used to find the electric field at x is 

E(x) = E (0 ) + — <N. - N ) , (2.7) 

e o 1 e 

where N, and N are the number of ion and electron sheets in the 
l e 

interval (0,x), and cj is the (positive) surface charge density of 
the sheet (see Fig. 2.2). Note that the electric field is discontinuous 
at the sheet position, and is not defined by Eq. (2.7). The electric 
field at the sheet position is taken at the middle of the jump shown by 
a cross in Fig. 2.2. Using the electric field thus computed, the Newton 
equation of motion for a sheet is integrated to find its position and 
velocity after a small time-step. From the new distribution of the 
charge sheets, the electric field can again be calculated, and 
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FIG. 2.1. Electric field generated by an isolated 
charge sheet. <j is the surface charge density. 



FIG. 2.2. Electric field calculated by Eq. (2.1). 
E(0) s 0 is assumed. 
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substituted in the Newton equation. Starting with a given initial 

distribution of charge sheets, the behavior of the system can thus be 

followed in time. When a particle leaves one end of the system, it is 

reintroduced at the other end, with the same velocity that the particle 

had upon leaving the system. This imposes a periodic boundary condition. 

Another condition sometimes used is that of a reflecting boundary: a 

particle rebounds elastically at the ends of the system. This boundary 

condition can be used for study of standing waves. 

In simulations of phenomena occurring on time scales for which 

only electron motions need be considered, the positive ions behave as 

an immobile neutralizing background. There are two ways of treating the 

ions in this case. One is to place ion sheets equidistant ly in space, 

5 

The other is to spread the positive charge uniformly. In the latter 
case, the electric field resulting from equally spaced electron 
sheets is as shown in Fig. 2.3, 


E 



FIG. 2,3. Electric field resulting from equally spaced 
electron sheets and uniformly distributed positive 
ion charge. 
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Simple as it is in principle, the sheet model described above is 

capable of simulating a wide variety of electrostatic, one-dimensional 

problems. However, since the number of particles that can be handled 

on modern computers is still many orders of magnitude smaller than the 

numbers of charged particles in experimental plasmas, the fluctuations 

appearing in the averaged quantities, such as the electric field and 

mean particle velocity, are much larger in a simulation model than in 

the plasma simulated. The fluctuation amplitudes can be reduced by 

increasing the number of simulation particles, N , but they decrease 
- 1/2 

only as N . Consequently, the energy of the fluctuations in a 

simulation plasma is typically many orders of magnitude higher than 
that in the plasma simulated. Since the enhanced fluctuations may 
cause nonphysical or exaggerated effects, considerable effort has been 
made to reduce the fluctuations without increasing the number of 
particles. This topic will be discussed in Section 2.3 in connection 
with various more sophisticated models developed from the simple charge 
sheet model just described . 

In two-dimensional problems, infinitely long charged rods are 
postulated, instead of charge sheets, to allow simulation in two 
dimensions. In three dimensions, the particles are represented by 
point charges. The electric fields created by isolated rod and point 
particles are shown in Fig. 2.4. 

In electromagnetic problems, for which the time-varying magnetic 
field is important, a model must be employed which properly solves the 
equation of charge particle motion, including a Lorentz force term, 
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dv . q . djx . 

XT" = — (E + v. x B) , = v 

d t m . ~ —1 — at 

1 


( 2 . 8 ) 


and the Maxwell equations [Eq. (2.2)], where subscript i denotes the 
particle species. The curi'ent and charge densities are obtained from 
knowledge of the positions and velocities of the particles. 



- r 



FIG. 2.4. Electric field, E , vs. radial distance, r. 

(a) Rod charge. q is charge per unit length. 

(b) Point charge. q is charge. 

P 

In what follows, we shall restrict ourselves to electrostatic 
phenomena, and concentrate on a one-dimensional model without external 
magnetic field. Therefore, we need only consider the one-dimensional 
equation of motion, 


dv . 


i 


dt 




and the Poisson equation, 


dE 

dx 


= 


(2.9) 


( 2 . 10 ) 
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Having considered the two approaches, it is an attractive possi- 
bility to develop a method which possesses the desirable features of 
both of them. Specifically, it should provide a very low fluctuation 
level, and retain particle discreteness effects, so that it will be 
possible to study both small amplitude waves and large amplitude non- 
linear phenomena. Such a 'hybrid* approach has been proposed by 
0 

Denavit. The main purpose of our work has been to develop this method, 
and to apply it to the study of nonlinear phenomena associated with 
large amplitude electron plasma waves. Section 3 is devoted to analyti- 
cal and numerical study of its characteristics. 

Various methods of solving the Vlasov equation are described in 
Section 2.2, Detailed descriptions of a variety of particle simulation 
models are given in Section 2.3. 

2.2 Solution of the Vlasov Equatio n 
2.2.1 Finite Difference Methods 

Since the Vlasov equation is a partial differential equation, 
probably the most straightforward approach is to approximate it by a 
finite difference equation, and to solve it numerically on a computer. 
This is simple, because it is a direct integration scheme, but lengthy, 
because a large number of grid points are necessary. 

The Vlasov equation is nonlinear since the electric field, E , is 
a function of the distribution function, f . It is a hyperbolic equation 
with variable coefficients. From the theory of numerical analysis, it 
is known that a finite difference method for a linear, hyperbolic, 
partial differential equation with constant coefficients is stable for 
appropriate choice of grid size, At and Ax , and convergent to the 


11 



exact solution in the limit -• 0 , fix - 0 . 


However, there is no 


guarantee of stability or convergence in the nonlinear, variable 


coefficient case. 

Consider the difference equation corresponding to Eq. (2.6) on a 

grid with x p = pax, t R = nAt, where p and n are integers, and 

f 11 _ f(x t ) A simple difference approximation to Eq. (2.6) is 
p P’ n 

obtained by the substitutions 


df 

at 



f - f 

5f P+1 Eli (2.11) 

ax 2 ax 


This yields 


n f /f “ f \ 

,n+l p+l + P-1 _ A1 [ P+- 1 Eli) (2.12) 

f p = 2 A x y 2 / 

which is known to be stable for jvAt/A x | < 1. The error in this 
scheme is of second order in At and Ax . It should be noted that 
the stability criterion is derived assuming that the coefficient, i.e. 
v is constant, so this must be considered as a local condition in 
velocity-space . 

A popular method of calculation, because it gives higher order 
accuracy, is the "leap-frog" scheme given by 




(2.13) 


which has the same stability condition jvAt/Ax| < 1 • The error in 
this scheme is of third order in At and A x . However, this has the 
disadvantage of being a three- time- level equation. A similar leap-frog 
scheme is used for advancing particles in the particle simulation models 


described in Section 2.3. 
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The nonlinear term, Egf/^v ( in the Vlasov equation can be treated 
in a similar way to the v$f/$x term. However, when the E&f/&v term 
is included, the difference equation may not be stable, even if the 
free- streaming part is approximated by a stable difference scheme. In 
studying two-stream instability for electrons streaming through mobile 
ions, Kellogg used a stable difference scheme for Eq, (2.6), but did not 
succeed in finding a stable scheme for Eq. (2.4), because of the effects 
of the nonlinear term, E^f/gv . 2 Various more sophisticated schemes 
have been pursued to minimize and suppress the effects of the numerical 
instability, including smoothing of the distribution function, and 
have been used successfully. 8,9 However, no general method of circum- 
venting this instability has been established; it depends very much on 
the nature of the problem under study. 

The difference equation for the Poisson equation [Eq. (2.5)] is 


E° = E n + 
P+ 1 P 


tX 


p.q 


1 q 

where Av is the velocity grid spacing, E* 1 = E(pAx, nAt), 

P 

n 

p q ~ ^ (pA x t qA v t n At ) , and the subscript i in f , denoting 
species, has been suppressed. For certain boundary conditions, 
convenient to solve 


(2.14) 


particle 
it is 



(2-15) 


where E = - d0/dx , A difference approximation to Eq. (2.15) is 


n 

^FH-1 



n 

+ *D- 


P-1 


(Ax) 2 


n 


(2.16) 
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The solutions of Eqs * (2.14) and (2.16) are straightforward, and cause 
no numerical instability. 

It is worth noting that the difference schemes mentioned above are 
"explicit" schemes, i.e., all the variables at the n-th time-step are 
expressed in terms of the variables at the (n-l)th time- step, or 
earlier. In an "implicit" scheme, a variable at the n-th time-step is 
expressed in terms of other variables at the same time-step and earlier. 
The use of the implicit scheme may make the difference equations stable, 
but since a set of simultaneous linear equations has to be solved at 
every time-step it may be prohibitively expensive. 

2.2.2 Transform Methods 

Transform methods constitute an obvious alternative to the 
finite difference method. Transformation of the variables to a different 
set offers the possibility of avoiding some of the difficulties 
encountered in the finite difference method. Expansion in terms of 
orthonormal eigenfunctions is a well-known mathematical technique, 

Fourier transformation being one of the simplest and most extensively 
used . 

Among the advantages of using the transform method are the elimina- 
tion of the partial derivatives in x and v , and of integration in 
v . The resulting algebraic operations are much simpler to deal with 
numerically than differentiations and integrations. Another advantage 
is a relief from numerical instability. This depends on the choice of 
transformation, and there is a possibility that a different type of 
numerical instability is introduced. The transform method reduces 
Eqs. (2.4) and (2.5) to a set of ordinary differential equations in 
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time, or to a simpler set of partial differential equations, depending 
on the transformation employed. The ordinary differential equations 
can be solved much more easily than the partial differential equations. 
Even when one of the partial dif f erentiations is left in the equations, 
they are easier to handle from a numerical stability point of view than 
the original equations. 

In the following, three methods will be described. They are the 
Fourier- Fourier method, the Fourier-Hermite method, and the Power 
Series Expansion method. For simplicity, we shall consider a one- 
dimensional, electrostatic, homogeneous plasma, without externally 
applied fields, and assume that the ions form an immobile, neutralizing, 
background charge. For the rest of the present subsection, we shall 
use dimensionless units for length, velocity, time, etc. Equations (2.4) 
and (2.5) can then be written as 


d£ 

at 


+ 






(2.17) 


where x is in units of the Debye length, , v is in units of the 
thermal velocity, , t is in units of the reciprocal of the plasma 
frequency, l/oj p , and E and f are also dimensionless variables. 

Fourier-Fourier Expansion** The distribution function, f , is 
expanded in a Fourier series in x , and Fourier transformation in v 
is carried out according to the following relations 
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f(x,v,t) = ~ Y: exp(jnk n x) / dy exp(-ivy) F n (y, t) , 


n=:“C» 


(2.18) 


■/ 


F n (y,t) = I dv exp (ivy) i / dx exp(-ink^x) f(x,v,t) f 

0 


where k Q = 2tt/L , and L is the length of the system. The electric 
field, E , is expanded in a Fourier series in x , 

CO 

E(x,t) = V E (t) exp(ink x) , 

/ j n 0 

n=:-oo 

(2.19) 

L 

V*>-E / E(x, t ) exp(-ink o x)dx . 

0 

By use of Eqs. (2,18) and (2.19), we obtain from Eq. (2,17) 


I? v*’*) + ^oh v^‘> - k - 0 - ° • 

(2.20) 

- ink Q E n (t) = F n (0,t) , 
where C n (y,t) is defined b 3 ' 

a> 

V y,t) = El v 0 ^ • «- 2i > 

m=-'co 

Equations (2.20) and (2,21) constitute the system of equations 
which is to be solved on the computer. Note that the integral over the 
distribution function in the Poisson equation [Eq. (2,17)] has disappeared. 
The electric field component E^ is not determined uniquely by these 
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equations, but is determined by the boundary conditions. For a periodic 
system without external excitation, we may put E^ = 0 . 

Equation (2.20) may be solved by integration along characteristics 
given by 

t - t/ = -jji- <y - y') . (2.22) 

0 

These are straight lines in the (t,y) plane, as shown in Fig. 2.5. 



FIG. 2.5. Characteristics of Eq . (2.20) in the (t,y) plane. 


The simplest numerical scheme is given by 

F n (y,t) = F n (y - nk Q At, t - At) 

At 

+ ^- (y - nk Q At) C n (y - nk^t, t - At) . (2.23) 
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For an improved approximation, we may use the iterative formula, 
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F (S+1) (y,t) = F(y - nk Q At, t - At) 

+ [(y " »At> C Q (y - nk Q At, t - At) + yC^ s) (y,t)] , (2.24) 

where superscript s denotes the number of iterations carried out # 
Equation (2.24) is obtained by replacing the last term in Eq. (2.23) 
by the average of the values of at the grid points (t,y) and 

(t - At, y - nk^t). Note that values of y in Eqs. (2.23) and (2.24) 
are chosen to fall at the grid points. Values of y - nk^t do not 
fall at the grid points unless £y and £t are chosen such that 
Ay = k 0 A t * Therefore, some interpolation scheme must generally be 
used. 

Since it is not feasible to include an infinite number of terms 

in Eq. (2,18) in the calculations, nor to cover the whole y-space, the 

Fourier expansion is truncated at n and a cut-off for F (y,t) is 

m n 

introduced at y = ± y^ . This is equivalent to a smoothing of the 
distribution function, f(x,v,t), expressed by 


f(x,v,t) = 


/ 


w (v')dv' 
v 


/ 

-00 


w (x')f(x + 

X 


v + v 


t )dx ' (2.25) 


where f(x,v,t) is the smoothed distribution function, and w (x) and 

x 

w^(v) are weighting functions given by 

n 

m 

w x (x) =1+2 ^ cos nk Q x , w v < v > - sin y m v • (2.26) 

n=l 
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The choice of the cut-off values, n and y must be made so that the 

mm 

half-widths of the weighting functions are small compared with charac- 
teristic lengths and velocities in the plasma. The truncations do not 
cause numerical instability.^ 

3 10 12 13 

Fourier-Hermi te Expansion ? ? * Instead of the Fourier trans- 

form in velocity used in the previous method, a Gram-Charlier series 
in velocity may be introduced, 

00 00 

f(x,v,t) = ^ ^exp(ink Q x) exp 
n=-co rrt=0 


0 


h(v}Z (t) 
m mn 


(2,27) 


where the h (v) are the orthonormal set of Hermite polynomials defined 
m 

by 


h (v) 
m 


, , v m r 2 y . 
(-1) exp(v /2) 

r A /2 -. 1/2 

[ (2n) ! m! ] ' 



(2,28) 


Again = 2 tt/L , and L is the length of the system. Electric field 
is expanded in a Fourier series according to Eq. (2,19). Substitution 
of Eqs . (2.19) and (2.27) in Eq. (2.17), and use of the recursion 
relations and orthogonality properties of Hermite polynomials and of 
Fourier series, yields 


dt 


Z (t) + ink 
mn 0 



Z 


m-1 , n 


(k) + 


(ra + 1) 


1/2 


Z 


rh* 1 , n 



m 


1/2 


E E (t)Z , (t) = 0 , 

n-q m-l,q ’ 


q— -co 


(m = 1,2,3, ; n = 0, ±1, ±2, ...) , 


(2.29) 
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where we have put E^(t) = 0 . The condition for reality of f(x,v,t) 
requires that 


Z 

ran 


(t) 



(2.32) 


Combination of Eqs. (2 . 29)- (2 . 31 ) provides a set of nonlinear, 
ordinary differential equations for the Z^Ct) . These may be solved 

step by step, typically by the Runge-Kutta method, or some improved 

. . 14 

technique. 

As in the Fourier- Fourier method, the infinite series expressed by 

Eqs . (2.19) and (2.27) have to be truncated. However, unlike the Fourier- 

Fourier method, the truncation of the infinite series [Eqs. (2.27)] in 

velocity-space causes serious difficulties. Truncation of the Fourier 

series can be justified provided that the perturbation amplitude is 

small. However, suppose that the Gram-Charlier series is truncated at 

m = M . From inspection of Eq. (2.29), it is clear that the expression 

for dZ /dt depends on Z _ (t) . The truncation consequently 

mn' mfl,n 

introduces serious inaccuracy, and therefore instability, when Z (t) 

x , n 

l/2 

becomes large. Since large by the time t =s M ' /(nk^) , 

the results of the computation are valid only for shorter times. This 
imposes a severe limitation on the usefulness of this method. A few 

attempts have been made to remove it, and the reader is referred to the 

^ 3,10, 13 

relevant literature. 
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Power Series Expansion ; Consider an expansion of F n (y,t) in 


Eq. (2,18) in powers of y , given by 


• < 2 - 33 > 

m=0 

where the g^ are arbitrary. The exponential factor improves the con- 
vergence of the series if it is truncated in the y-direction. Substi- 
tuting Eq. (2.33) into Eq. (2.20), and collecting equal powers of y , 
gives 


— a - nk 
dt rnn 


‘ g ra-l g ra+l 

„ a , - (m + 1) 

°L s » ”- 1 '“ s « 


m+1, n 


m-1 

Vm 


a A a - = 0 , 
Oq m-l,n-q 


(2.34) 


q=-co 


(n = 0, ± 1 , ± 2, m = 0, 1 , 2, ...) . 


This equation can be solved in a similar way to Eq. (2.29) obtained in 

the Fourier-Bermite method. It should be noted that the coefficients 

a may be shown to be equal to the Z in the Fourier-Hermite method 
mn M mn 

15 

except for a complex factor. 

In the numerical integration, Joyce et al . have encountered the 

difficulty associated with the truncation of the infinite series in 

15 

velocity space that was found in the Fourier-Hermite method. However, 
they found that numerical solution of the linearized version of Eq . (2.34) 
gave a very regular pattern for the amplitudes for large m , as 

shown in Fig. 2.6. The coefficients seem to be sampled points of a 
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continuous function in m . This would lead to a guess that the 

(m +1) coefficient may be extrapolated by polynomials, and that it 
max 

may be possible to close the system expressed by Eq. (2.34). This 
method turned out to work well for the linear Vlasov equation, but 
Joyce et al. found that for nonlinear cases its success depends on the 
problems treated. Thus, although there is improvement over the 
Fourier-Hermite method, in the sense that it stabilizes the integration 
scheme, the method is still not satisfactory for simulations carried to 
long times. 


a m xl0 3 



FIG. 2.6. An example of amplitude a vs. m . 

mn 

(Adapted from Fig. 1 of Ref. 15.). 
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2.3 Particle Simulation Model 


A number of authors have developed schemes which reduce the high 
level of fluctuations characteristic of the electrostatic sheet model, 
without affecting the collective behavior of the plasma. 16-26 The 
reduction of fluctuations in these schemes is achieved by smoothing the 
interaction forces between the charge sheets. The smoothing reduces 
the short wavelength fluctuations arising from discontinuities in the 
electric field generated by the sheets, but leaves collective phenomena 
associated with long wavelengths unaffected. Since the smoothing may 
modify the physics of the plasma, and perhaps introduce spurious effects 
it is important to examine very closely the detailed behavior of the 
model. In this section, we shall describe various electrostatic plasma 
models incorporating smoothing procedures, and discuss some of the 
physical consequences of the smoothing relevant to the hybrid approach 
to be treated in Section 3. All of these schemes can be used for multi- 
dimensional simulations, with or without an externally applied magnetic 

field, but the discussion here will concentrate on the one-dimensional 
case without magnetic field. 

Zero-Size-Particle and Nearest-Grid-Point Method 

In this method, instead of treating the space variable, x , 
as a continuous variable when computing the electric field, the system 
is divided into an arbitrary number of cells, and the electric field is 
computed only at grid points corresponding to the centers of the cells. 
The charge distribution from which the electric field is calculated, 
through the Poisson equation, is obtained by accumulating the charge 
of all the sheets in each cell at the center, regardless of their 
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positions within the cell. In advancing in time, all of the particles 

in one cell are accelerated by the electric field defined at the grid 

point, regardless of their positions. Note, however, that the particle 

position is not discretized, but is a continuous variable. This scheme 

was used by Burger et jal. in a study of randomization mechanisms in a 

1 6 

collisionless plasma. Hockney applied it in a two-dimensional simu- 
lation of anomalous plasma diffusion across a static magnetic field, 

17 18 

using a charge rod model. * 

In one dimension, the electric field calculation is quite simple. 
The charge distribution on the grid points is obtained by counting the 
number of electron and ion sheets separately in individual cells, 
multiplying them by the surface charge density, and subtracting the 
results. The electric field at a given grid point is then taken at the 
middle of the jump in electric field occurring at that grid point. When 
ions can be assumed to be immobile, equal numbers of ion sheets are 
assigned to each cell, and only the number of electron sheets needs 
to be counted , 

Advancing particles in time through the Newton equation may be done 
by the leap-frog scheme. In one dimension this may be written, for an 
electron sheet, as 


V n+l/2 “ V n-l/2 



At , 


X 

n+1 


x 

n 


+ v 


n+1/2 


At 


( 2 . 35 ) 


where v , and v 

n+1/2 n- 

half time-step ahead of, and behind, the n-th time-step, respectively, 


x and E are the sheet position and electric field at the n-th 
n n 


1/2 


denote the velocities of the sheet at a 
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time step, and e are electron mass and charge (magnitude), and 

is the time- step . This scheme is widely used because of its high order 
accuracy [0 { (At) 3 }] . 

The ZSP-NGP scheme not only reduces fluctuations, but also 
makes the computation much faster than for the simple sheet or rod model. 
Furthermore, the gain in computational speed is made by using integer, 
rather than floating-point, arithmetic on the computer. In most 
computers the integer addition and subtraction are faster than the 
f loating-point addition and subtraction. Since the main part of the 
algorithm can be written in terms of additions and subtractions, computing 
time can be saved by normalizing and expressing the particle positions 
and velocities, and the charge and electric field, in terms of integers. 
The reader is referred to the relevant papers for more details . ^ 9 

The reduction of fluctuations can be understood as follows; the 
interaction force between a charge sheet located at a grid point, and 
another sheet separated by distance, x , will be as shown in Fig. 2.7(a). 
The method eliminates forces acting in close encounters, whereas in the 
simple sheet model there is no zero-force region. This helps to reduce 
collisional effects, and fluctuations at short wavelengths, but does not 
affect the collective behavior due to long range forces. 

The force between two charge rods is illustrated in Fig. 2.7(b). 

It is a function varying in steps. At large separations, it is a good 
approximation to the force between two line charges. At close separa- 
tions, however, the force vanishes, as in the one-dimensional case, 
whereas the force between two line charges would tend to infinity. 
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FIG. 2.7(a). Interaction force in ZSP-NGP model, between a 
charge sheet at a grid point, and another sheet at distance 
x apart. (Adapted from Fig. 5 of Ref. 29.). 
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FIG. 2.7(b). Interaction force in ZSP-NGP model, between two 
charge rods distance r apart, r is in units of cell 
size h . (Adapted from Fig. 1 of Ref. 17.). 


Although the ZSP-NGP method has been used successfully, 1 ^ ^ the 
stepped force law shown in Fig. 2.7 may not be accurate enough, and the 
zero-force region may smooth out small amplitude oscillations, even for 
long wavelengths. The discontinuity in the electric field still produces 
undesirable enhanced fluctuations. One way to relieve these difficulties 
is to increase the number of particles and to use a finer grid. In 
practice, however, it is not always feasible to increase the number of 
particles and grid points to a satisfactory level. The finite-size 
particle models which we shall discuss next have been developed to over- 
come these difficulties, particularly in the simulation of high density 
plasmas. We shall describe the Cloud-in-Cell method, Particle-in-Cell 
method, Gaussian Cloud method, Multipole Expansion method, and Lewis’ 
method . 

2 . 3.2 Cloud-in-Cell and Particle-in-Cell Methods 

It is natural to introduce the idea of charged particles of 
finite dimensions, for the purpose of improving upon the ZPS-NGP method. 

In the Cloud-in-Cell method, the discontinuities are mitigated by 
considering clouds of uniformly distributed charge which are assumed to 
be tenuous, and to be able to pass through one another. 19 The center 
of the charge cloud is taken to be the particle coordinate, and a 
spatial grid is used for computing field quantities. Key features 
of the method are the prescription for determining the amount of charge 
to be assigned to the grid points (charge-sharing), and the force on a 
cloud from field quantities on the grid points (force-sharing). The 
charge is distributed to neighboring cell centers (grid points) in 
proportion to the areas the cloud occupies in the cells. The charge 
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density on the grid points thus accumulated is used to compute the 
electric field. The force on a cloud is calculated by averaging the 
contributions from the neighboring grid points with weightings propor- 
tional to the areas in the corresponding cells occupied by the cloud. 
For example, in the one-dimensional case shown in Fig. 2.8(a), the 



(b) TWO DIMENSIONS 


FIG. 2.8. Finite-size particle located in a grid. Shading shows 
assignment of charge density to grid points in CIC model. 

(a) Ax is cell size. H is particle size. 

(b) ^ and Ay are cell sizes. A is area of a particle. 


charge assigned to each grid point is given by 



“i - *i + i - H ■ < 2 - 36 > 
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where p is the charge density of the cloud, H denotes the size of 


the cloud, and a^ represents the portion of the cloud occupying the 
i-th cell. The extension to the two-dimensional case is shown in 
Fig. 2.8(b). The charge assigned to each grid point is given by 


■t., ■ 4 ¥) ■ p 1+1>J ■ 4 ^) ■ 


p i+l,j+l P c 




P i,j+1 P c 


M ■ 


(2.37) 


a. a, "{”3 “hft = A 

i,j i+1, j T i+l.J+l i| j+1 


where A denotes the area of the cloud, and a. represents the 

1 > J 

portion of the cloud occupying (i,j) cell. 

The cloud size may be either larger or smaller than, or equal to, 
the cell size. It is chosen so that the particle dynamics can be 
adequately represented for the problem of interest. This will be clear 
from Fig. 2.9, which shows the density assigned to grid point i for 
clouds of various sizes moving in the x-direction. As the cloud 
size increases beyond the cell size, the resolution decreases, because 
the cloud density is held constant over a distance larger than the 
shortest resolvable wavelength, i.e. the cell width. 

The force on the cloud is given by, 

a. 


' - '.Sir'i 


Z 3 p< 


Z.Yi 


(one dimension) 


(two dimensions) 


( 2 . 38 ) 


ij 


id 


29 




FIG. 2.9. Sketch of charge density assigned to grid point 
i as a finite-size particle moves in x-direction. The 
size varies from 0 to 2£x . The horizontal axis repre- 
sents the position of the center of the cloud. (Adapted 
from Fig. 3 of Ref. 19.). 


where E^ (or IS^) electric field at the grid points (see 

Fig. 2.8). This force-sharing scheme produces no self-force if it is used 
together with the charge- sharing scheme just described. If this charge- 
sharing scheme were not followed, the force averaging could lead to 
incorrect acceleration of clouds. 19,29 

The interaction force between a cloud fixed at a grid point, and 
another cloud is plotted as a function of separation in Fig. 2.10, It 
will be noted that the interaction is much smoother than in the ZSP-NGP 
method, and that the zero-force region is eliminated. The CIC scheme 
consequently produces less collisional effects than the simple sheet 
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FIG. 2.10(a). Interaction force in CIC model, between a 
cloud at a grid point and another cloud at distance x 
apart. [H = fix] . (Adapted from Fig. 6 of Ref. 29.). 
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FIG. 2.10(b). Interaction force in CIC model, between a 

positive charge cloud at a grid point and a negative charge 
cloud, at distance x apart but with the same y 
coordinate. [H = Ax]. (Adapted from Fig. 4 of Ref. 19.). 
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model, and is better able to simulate small amplitude oscillations. It 


has been observed by Birdsall and Fuss that the CIO method, combined 

with the leap-frog algorithm, greatly improves the energy conservation 

19 

over the ZSP-NGP method for the same number of simulation particles. 

In the Particle-in-Cell method, proposed by Morse and Nielson, the 

charge is assigned to each cell by distributing the charge of a particle 

20 21 

between the nearest two cells according to a linear interpolation. , 

The electric field is then calculated from the resulting charge distri- 
bution on the grid points, in the same way as in the ZSP-NGP and CIO 
methods. The force on a particle is obtained by calculating the electric 
field at the particle position by linear interpolation between the 
nearest two cell centers. It is easily seen that this scheme is a 
special case of the CIC scheme, with the particle size equal to the 
cell size. 


2.3.3 Gaussian Cloud Method 

Instead of a uniform charge dis tribution , in this method a 
cloud is supposed to have a Gaussian charge distribution, 


Pi<x) 


Pc 


exp 



2a 


(2.39) 


where x i is the position of the center of the cloud i , and a 
represents the one-dimensional size of the cloud. 

To compute the electric field from a given distribution of charge, 
it is convenient to work in Fourier- transformed space rather than in 
coordinate- space, i.e. in k-space rather than x-space. The Poisson 
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equation then simplifies to an algebraic form, and is handled more 

easily. The Fast-Fourier-Transform technique can be used to speed up 
30 

the computation, 

A model proposed by Dawson et al. makes use of the foregoing ideas 
22 

in the following way. The electric field generated by a Gaussian 
cloud is obtained by the Fourier transform, 

E ( x ) = /■<*, exp(ikx)dk , (2.40) 

where E(k) is given by 


4np 

ikE(k) = -S- exp 

(2tt) /2 


,2 2 

k a 


+ ikx . 


(2.41) 


The force on particle i due to particle j is then given by 


F iJ =/V X) Pi (x)dX 

= 2p c ex PC ik < x i “ x j) - k 2 a 2 ]dk . (2.42) 

The integral over k-space is replaced by the finite sum 

k 

max 

F ij = 2p cS E ex PC - k 2 a 2 ) sin k(x. - x . ) , (2.43) 

k . 
min 

where k = 2 n n/L , n is an integer, and L is the length of the 
system. To advance particle i in time, the total force on it is cal- 
culated by summing the contributions from all other particles, and is 

assumed to be constant during a time-step. The leap-frog scheme may be 

23 

used in this model. As will be shown in Section 2.3.6, the 
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fluctuations with short wavelengths can be reduced substantially below 
those in the simple sheet model as the value of a increases. 


2.3.4 Multipole Expansion Method 

The multipole expansion scheme was developed by Kruer and 
24 

Dawson. It was used by Denavit and Kruer, and compared with numerical 
solutions of the Vlasov equation for two-stream instability and sideband 
instability.^ Close agreement was obtained between the results of 
the two methods. Consider a charge cloud with an arbitrary charge 
distribution and a spatial grid. The location of the center of the 
particle is defined by 

x i = x Gi + 6x^ , (2.44) 

where x , is the nearest grid point location, and 6x. is the particle 
Gi l 

coordinate relative to x . The Fourier transform of the cloud charge 

Gi 

density, p^(x), of particle i is 


P .(k) 



exp(-ikx)dx 


= S(k) exp(-ikx„. ) exp(-ik6x.) , (2.45) 

Gi 1 

where S(k) is the form factor, which is the Fourier transform of the 
shape factor given by 


S(x - x . ) = p . (x) . 

l F i 

Assuming kfix^ « 1 , we may expand p^(x) as follows 


(2.46) 


D , (k) = S (k) exp(-ikx . ) (1 - ikfix. +...). (2.47) 

K l Gi l 
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This is equivalent to replacing the particle centered somewhere in the 
cell by a particle centered at the nearest grid point, plus a dipole 
there and higher-order multipole terms. Summing over a collection of 
particles, and inverse Fourier-transf orming the Poisson equation, yields 
the electric field at the grid point. 

To compute the force on a particle, given by 



Physically, this means that the force on a particle may be represented 
by its monopole times the electric field, plus the dipole moment times 
the derivative of the field, plus higher-order contributions. 

The scheme just outlined is valid only if k fix. « 1 , In order 

max 1 

for this to be satisfied, the particle size, H , has to be larger than 
the cell size, &x , because fix i £ A*/ 2 and k m ax H ~ 1 • In P rac tice, 
the expansion is truncated at the dipole term, and a particle size of 
several cell-lengths is necessary, i.e., if the particle size is 
comparable to the cell size the simulated behavior of short wavelength 
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perturbation is likely to be inaccurate. When the short wavelength 
components are not important, it may be adequate to use this scheme 
even for H ~ £x . 

2.3.5 Lewis 1 Method 

In the simulation of plasma phenomena, it is important to 

check whether the model conserves energy. If there is no external 

excitation, the total energy of the system should remain substantially 

constant during the simulation run. Otherwise, it is to be suspected 

that nonphysical effects are being introduced. There is no guarantee 

that the models described so far are energy-conserving. Lewis has 

presented a general method for deriving numerical approximation schemes 

25,26 

which guarantee energy conservation, 7 It is based on a Lagrangian 
formulation via Hamilton’s variational principle. The CIC method is 
closely related to a special case of the algorithms based on this 
formulation. The Lewis method serves to establish both the theoretical 
basis for the use of various models in plasma simulation, and to pro- 
vide an alternative viewpoint on the problems involved. We shall 
outline the method and show how the energy-conserving version of the 
CIC scheme is derived by Lewis. 

For simplicity, we consider a one-dimensional, electrostatic 
plasma with no externally applied fields and neutralizing immobile ions. 
The Lagrangian for this system is 


/ ( m X ) 1 f 2 

dx 'dv'f (x', v', 0) | — + e0(X, t) j ^ J dxE~(x, t) , (2.51) 


where X[= X(x',v',t)] describes particle trajectories as functions of 
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the time, t , initial position, x' f and velocity, v' , the dot denotes 
differentiation with respect to t , f(x,v,0) is the initial distribution 
function, 0 is the electrostatic potential, and E = -30/ , We 
approximate 0 and X by finite series 




N 


0(x f t) = ^ a n (t)§ n (x) 


n=l 


X(x',v',t) = 7 ^(t)X^(x', v') , 

1=1 

(2.52) 


where $^(x) and X^(x',v') are linearly independent basis functions, 

and $ (x) must satisfy the boundary conditions. Substitution of 
n 

Eq. (2.52) into Eq. (2.51) provides a new Lagrangian which is a 
function of the generalized coordinates, a n , , and the generalized 

velocities, y . Applying Hamilton's principle, 




6 / X dt = o , 

t. 


(2.53) 


for arbitrary variation of qi^ and 
equations , 




we obtain the Euler-Lagrange 


d_ 

dt 


(«;)- 


ax 

& 7 , 


= 0 , 


ax 


= 0 . 


(2.54) 


These are the equations which describe energy- conserving numerical 
approximation schemes. It remains to specify functions and X^ . 

Let the initial distribution function be a sum of the Dirac delta- 
functions , 
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(2.55) 


f(x,v,0) 6(x - x ^ ) s(v - v i > , 


where x^ and v are the initial position and velocity of the i-th 

particle. Choosing a set of X -functions as 

/ £/ 


X^x'.v') - 


U' , x t . v' - v t ) . 


(otherwise) , 


(2.56) 


implies that the initial conditions on the y (t) are 

■v 


V 0) = \ ' 


V 0) = v -t • 


(2.57) 


Clearly, y^(t) is the position of the l~th particle at time t . To 

obtain an energy-conserving version of the C1C scheme, it is convenient 

to take $ (x) to be 
n 


r-[x - (n-l)h] [(n-l)h < x < nh] , 


h 


$ Q (x) = ^ -i-[(n+l)h - x] [nh < x < (n+l)h], (2.58) 


h 

0 


[otherwise] , 


where h = l/(N.j+ 1), and L is the system length. With this set of 
functions, Q^Ct ) is the value of potential at x = nh . 

Substituting Eqs . (2.56) and (2.58) into the Euler-Lagrange equations 
[Eq. (2.54)] gives 


m 


h * . 


(2.59) 
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a ,(t) - 2a (t) + a n ,(t) 
n+ 1 n n~ 1 


Po 


€ o h 


£ ' (2 - 60) 


where $^(x) = (d/dx )$^ (x) , and is the background charge density. 

Equation (2*59) is the equation of motion for the £-th particle and 
Eq . (2.60) is the central difference approximation for the Poisson 
equation. It is seen that Eq . (2.60) is identical to the approximation 
for the Poisson equation in the CIC method. Equation (2.59) differs 
from the equation of motion in the CIC method only in that a piecewise 
linear function is used in place of (t)] . The CIC method con- 

l -v 

sequent ly does not conserve energy. Equations (2,59) and (2.60) are 
the energy-conserving version of the CIC scheme. It should be noted 
that the energy conserving feature of this method is realized only in 
the limit of a vanishingly short time-step. In practice, the time-step 
is small but finite, so the energy conservation is not exact. 


2.3.6 Finite-Size Particle and Spatial Grid Effects 

In this section, we shall consider the effects of intro- 
ducing finite-size particles, and a spatial grid into simulation models. 

31 29 

We shall follow the analyses given by Langdon and Birdsall and 
consider, first, how longitudinal plasma oscillations are affected, and 
second, how the fluctuations are reduced. 

In general, finite particle size enters the analysis via the form 
factor, S (k) , i.e. the Fourier transform of the shape factor of the 
charged particle distribution, S(x) , 

S(k) = /dxS(x) exp(-ikx), (2.61) 
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Using the form factor S(k) , the charge density, current density, and 


force in the finite-size particle model can be written as, 


p c (k,t) = S(k) p (k,t) , 


J <k,t) = S(k)J (k, t) , 
c p 


F (k,t) = S(-k)F (k,t) , 
c p 


( 2 . 62 ) 


where subscript p in p , J , and F is affixed to emphasize 
that they are for a system of point particles. These relations suggest 
that the theory of finite-size particles may be obtained by multiplying 
the charge in the point particle theory by S(k) . This is indeed so 
for the longitudinal plasma permittivity: for a one-dimensional plasma 

with finite-size cloud, we have 


2 % r a V av 

= 1 * S <k) / — n- 


(2*63) 


where Landau’s prescription applies regarding analytic continuation. 

For a Maxwellian velocity distribution, and CIC-type finite-size 
particles, Eq . (2.63) becomes 


€ p ( k ,U)) 
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S 2 (k) 
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2k X D 




S(k) 


sin (kH/2) 
(kH/2) 


(2.64) 


where H is the size of the cloud, and Z ' denotes the first derivative 
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of the plasma dispersion function. For a Gaussian cloud, such as is 
used in the Gaussian Cloud method, the form factor becomes 

S(k) = exp(-k 2 a 2 /2) . (2.65) 


40 



Note that in the limit of vanishing size, S( k) tends to unity and the 
point particle theory, i.e. the simple sheet model, is retrieved. 

We have obtained solutions of the dispersion relation [ e (k,tn) = 0] 
for point particles, and for finite-size particles. They are plotted 
in Fig. 2.11. It is clear that the long wavelength oscillations are 
little affected, while the short wavelength oscillations are strongly 
modified . 

When spatial grid effects are taken into account, the permittivity 
31 

is given by 
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S e (k) = S(k) 


sin (kAx/2) 
(k^x/2) 
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sin (kAx/2) 
(kAx/2) 
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k 


sin k£x 
k A x ’ 



( 2 . 66 ) 


where l is integer. 

This dispersion relation has been studied in detail by Langdon, 31 
who has shown that small values of \^/ fcx lead to instability. This 
is due to the coupling between waves with different wavelengths through 
the spatial grid. He has also shown that when X^/^x s 1 , the insta- 
bility is negligible, and the = 0 term in the summation in Eq. (2.66) 
is dominant. The result of our computation for dispersion characteristics 
with both finite-size effects and spatial grid effects taken in account 
is shown in Fig. 2.11, Only the L = 0 term has been retained in 
Eq. (2.66). It will be seen that the spatial grid exerts an additional 
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FIG. 2.11. Solutions of dispersion relation. 
[Eq. (2.66)]. 
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smoothing effect on the dispersion characteristics, modifying the 
behavior of short wavelength oscillations in a similar way to finite- 
size effects, leaving the collective behavior at long wavelengths 
almost unchanged. 

Finally, we turn our attention to the thermal fluctuation spectrum. 
The fluctuation spectrum of the electric field in k-space for a 
point particle plasma is 



(2.67) 


where is the Boltzmann constant, and T is the temperature of the 

plasma. Neglecting the coupling between perturbations due to the grid, 

the fluctuation spectrum for a cloud plasma in a gridded system is 
33 

given by 



k B T / k 2 / 1 

2 (k 2 ) ( K /k + K 2 X^(k)) 


( 2 . 68 ) 


where k , K , and are given in Eq. (2.66). In Fig. 2.12, the 

spectrum computed from Eq. (2.68) is plotted for the CIC finite-size 
particle model, and for the Gaussian Cloud model in a gridless system 
(£x -* 0). The thermal fluctuation spectrum for the point particle 
plasma [Eq . (2.67)] is also shown for comparison. Note that the 
comparison in Fig. 2,12 is for the same number of particles. The 
reduction of fluctuations in the short wavelength part of the spectrum 
is apparent. An interesting feature of the plots is that the CIC 
spectrum has zero-energy holes, corresponding to modes satisfying 

S(k) = 0 , while the Gaussian Cloud model does not. 

43 




0.1 I to 

kX[) 


FIG. 2.12. Fluctuation spectrum of electric field 
[Eq. (2.68)]. 



The modification due to spatial grid effects in the C1C model is 
also shown in Fig. 2.12. It will be seen that introducing a spatial 
grid reduces the short wavelength fluctuations further, and creates 
additional holes in the spectrum. 

The analytical results of Eqs. (2.66) and (2.68) were verified 

by Okuda in computer studies of the CIC model. 34 Okuda also performed 

a series of simulations to study the numerical instability due to a 
35 

spatial grid, and confirmed that the simulation is stable when the 
Debye length is comparable to, or larger than, the grid size, as 
predicted by Langdon. 31 

It is evident from the foregoing considerations that finite-size 
particle models are capable of simulating the collective behavior of 
plasma, with reduced fluctuation level, by choosing appropriate particle 
and grid sizes. However, the reduction in the fluctuations in these 
methods is not more than about 10 dB below the simple sheet model. It 
is consequently still very noisy compared with the practical plasmas 
being simulated. In Section 3, we shall describe a method which provides 
a fluctuation level several orders of magnitude lower than the finite- 
size particle models. 
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3. LOW-KOISE HYBRID APPROACH 


3 . 1 Comparison between the Vlasov and Particle Simulation Approaches 
Before describing the hybrid approach used in our work, it is 
convenient to compare here the two basic approaches to computer solution 
of plasma problems discussed in Section 2. We shall do so by considering 
how well the plasma dynamics can be represented, and how easy the com- 
puter solution is to accomplish. 

3.1.1 Plasma Dynamics 

The Vlasov approach brings into question the validity of the 

Vlasov equation itself: the Vlasov equation is a description of a 

plasma which is correct only to lowest order in the plasma parameter, 

3 -1 

(nXp) , where n is the number density, and X^ is the electronic 
36 

Debye length. As a consequence, it describes collective effects due 
to long-range Coulomb forces, but does not include particle discreteness 
effects such as particle-particle encounters. Since its solution con- 
tains no thermal fluctuations of macroscopic quantities, such as the 
electric field, charged particle density, etc., the behavior of very 
small amplitude waves can be simulated in a quantitatively accurate 
manner , 

In contrast, the particle simulation model, or ’particle code’ as 
it is often termed, incorporates the full dynamics of particles and waves, 
including discreteness effects. However, fluctuations are at an 
unrealistically high level, since the number of particles that a computer 
can handle is many orders of magnitude less than in typical plasmas. 

These large fluctuations can obscure the small amplitude oscillations 
described by linear theory. Hence, quantitative results are not easily 
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available in the linear regime, even though the qualitative behavior 
agrees with theory; the particle code is better suited to simulation in 
the nonlinear regime, where the ratio (wave field energy/particle thermal 
energy) > 0.01. In this respect, it should be noted that the expansion 
methods discussed in Section 2.2.2 may limit the amplitude of the 
perturbations that can be handled in order to assure the rapid conver- 
gence of the expanded series, it should also be noted that the expansion 
methods impose a restriction on the class of perturbations that can be 
treated. For example, in the Fourier-Hermite method of Section 2.2.2, 
the velocity distribution function is automatically an entire function 

3 

of velocity. Consequently , the problem must be such that the 
distribution function is describable in terms of an entire function of 
velocity; we could not consider the evolution of a delta- function 
distribution . 

3.1.2 Ease of Computation 

From the computational point of view, it is much more 
difficult to solve the Vlasov equation numerically than to solve the 
Newton equation. Numerical instability such as is encountered in the 
finite difference method, and the Fourier-Hermite method, is difficult 
to overcome. Solution of the Newton equation is not subject to this 
difficulty. In addition, the algorithm for solving the Vlasov equation 
is more complicated and lengthy than for the Newton equation. 

One point of practical importance is that the Vlasov approach 
enables different effects to be studied separately. For example, the 
linearized Vlasov equation can be solved to study linear behavior alone. 
The effects of nonlinearity can then be assessed by adding the nonlinear 


47 



terms. The effect of free- streaming particles can also be studied 
independently. This feature of the Vlasov approach helps materially 
in understanding the physics of plasmas . 

The particle code requires a very large computer memory to store 

information on the positions and velocities of the particles, and on 

the field variables, since it follows the motion of the particles 

step-by-step in time. Even with modern, large, high-speed computers, 

simulations in three dimensions may not be economically feasible, since 

the number of particles that is needed, to maintain the same level of 

3 

fluctuations as in one dimension, has to be increased as N 

In contrast, the Vlasov approach needs much smaller computer 
memory for the time-varying field quantities on the grid points, or the 
time-varying expansion coefficients, since the velocity distribution 
carries information on the charged particles in highly compact form. 

To ease the computational difficulties, some combination of the 
two approaches would clearly be desirable. This is provided by the 
hybrid approach to be described next. It is a particle simulation 
model in the sense that the motions of a large number of particles are 
followed in time. The particles do not keep their identities, however, 
and there are similarities to the Vlasov approach in that the value of 
the velocity distribution function is defined on a grid in phase— space, 
as in the finite difference methods described in Section 2.2.1. 

3.2 Hybrid Approach 

Our model is constructed using the CIC model described in 
Section 2.3.2. In addition to the usual spatial grid, we introduce a 
grid in velocity-space, and represent the particles by points in 
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(x v) phase-space as shown in Fig. 3.1. The phase-space is consequently 
covered with a rectangular grid of dimensions £x, . The velocity 

grid extends from to , where and v^ are chosen such 

that the numbers of particles with velocities in the intervals 
v < v^ or v > are negligible. 

In the work of Denavit,*^ the Lewis variational method 2 ^’ 26 was 
used to construct a model. In the Denavit model, finite-size particles 
are chosen to have a triangular spatial distribution, instead of the 
uniform charge distribution of the CIC model. The numerical scheme 
based on that model turned out to be more complicated, therefore more 
time-consuming, than our model. In addition, although his scheme is 
energy-conserving, it does not conserve momentum, whereas the CIC scheme 
is formulated in such a way as to conserve momentum. 2 '*' The non-conserva- 
tion of momentum indicates the existence of self- force, i.e., a particle 
is effected by the force due to the field that is created by itself, 
which is nonphysical. The energy-conserving feature of the Denavit 
scheme may not be very useful in practice, since energy conservation is 
exact only in the limit of a vanishingly small time-step. 

In creating a plasma with a Maxwellian velocity distribution, our 

model employs the following method. The particles are equally divided 

into a number of velocity groups. All the particles in one group are 

assumed to have the same velocity, v , and mass and charge are assigned 

to them in proportion to exp (-v 2 /2v 2 ) , where v is the electron 

t t 

thermal velocity. This is shown schematically in Fig. 3.1. Since the 
charge-to-mass ratio is the same for all of the particles, the accelera- 
tion is also the same. One of the advantages of this method of 
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HO. 3.1 . Phase-space covered with a rectangular grid, and a Maxwellian velocity 
distribution approximated by beams. 


generating a Maxwellian distribution by weighted particles is that 
improved resolution is provided in the tail of the velocity distribution, 
compared with a Maxwellian distribution with identical particles. 

The system is set up at time t = 0 using the quiet start technique 
to be described below, and proceeds as in a CIC model. After a certain 
number of time-steps, the distribution function is reconstructed at the 
grid points by periodic smoothing, and is interpreted as representing 
a distribution of new discrete particles. The motions of these particles 
are followed until the next reconstruction. 

With these procedures in mind, we shall describe the quiet start 
technique, and the periodic smoothing, which are essential parts of the 
hybrid approach. As proposed by Denavit, they are used to achieve a 
very low fluctuation level, and to allow the model to be applied to a 
wide range of linear and nonlinear problems. 

3.2.1 Quiet Start 

The quiet start technique, proposed by Byers, is a method 

of eliminating macroscopic fluctuations in a particle code at early 

37 

stages of evolution. This is done by placing the particles only on 
the equilibrium trajectories in phase-space at time t - 0 . As an 
example, consider an electrostatic, one-dimensional problem with a 
periodic boundary condition, in which the equilibrium particle trajec- 
tories are straight lines, v = const. At time t = 0 , the particles 
are loaded uniformly at the grid points in phase-space, as shown in 
Fig. 3.1. We shall assume that the particles are of finite extent, as 
in the CIC method of Section 2.3.2, with extension equal to the spatial 
grid size, Ax . Since the charges are distributed uniformly in 
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space, there is no macroscopic electric field at t = 0 . We easily 


see that at subsequent times the phase- space looks the same as at 
t = 0 , except that the centers of the finite-size particles are shifted 
in the x-direction by distances depending on their velocities . Thus, 
there continues to be no macroscopic electric field unless a perturbation 
is applied. 

In principle, the quiet start technique provides a noiseless 
system. In practice, round-off errors due to the finite number of 
digits representing the numbers in the computer introduce some fluctu- 
ations. However, this level is many orders of magnitude lower than that 
in the particle codes. Typically, we have obtained about 60 dB reduction. 
This is quiet enough to study linear plasma behavior, and compare its 
features with theoretical predictions. 

Although the quiet start technique works well at early times, it 
eventually causes wave growth; it ceases to be effective after a time, 

6 

Av). where k is the maximum wavenumber possible in the system, 
m m 

This breakdown occurs because the velocity distribution of the particles 

is being replaced by a set of discrete beams. Such a system is subject 

to streaming instability, to be described in Section 3.3, even if the 
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envelope of the beam density is Maxwellian. Periodic smoothing may 
by used to combat this instability. 

3.2.2 Periodic Smoothing 

Periodic smoothing constitutes a periodic averaging of the 

0 

distribution function in phase-space. It can be expressed by 



f(x,v,t) =JJ f (x',v',t)w x (x - x')w v (v - v ' )dx 'dv ' , (3.1) 
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where w x and w y are weighting functions for coordinate- and velocity 
space, f is the averaged distribution function, and the integration is 
over the whole of phase-space. The averaging process expressed by 
Eq. (3.1) causes diffusion of the distribution function. The weighting 
functions are chosen so that this diffusion quenches the streaming 
instability without introducing other undesirable effects. Specific 
forms of the weighting functions, and their derivation, are given in 
the Appendix. 

In particle models with a phase-space grid such as that shown in 
Fig. 3.1, the integral in Eq. (3,1) reduces to a sum over the collection 
of particles, and we want to find the averaged distribution function, 
f, at the phase-space grid points, if f(x',v',t) is taken to be the 
mass of a finite-size particle, the center of which is located at 
(x'jV 7 ) at time t , then Eq. (3.1) implies that the value of f at the 
(i-j) grid point, (x..,Vj), is obtained by distributing the mass of each 
particle among the neighboring grid points according to the weighting 
prescribed by w x and w y . This is a reconstruction of the distri- 
bution function from a given distribution of particles. In this model, 
the reconstructed distribution function is defined only at the grid 
points, and is interpreted as the new particle mass located at each 
grid point. Their motions are governed by the Newton equation, as in 
other particle codes. Reconstruction of the distribution function does 
not need to be done every time-step. It is simply done frequently 
enough to suppress the streaming instability, as described in Section 3.3 
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3.3 Streaming Instability and Recurrence Phenomenon 


In this section we shall discuss two phenomena associated with 
the use of the quiet start technique: streaming instability, and 

possibility of recurrence of the initial state. 

3.3.1 Streaming Instability 

Treating each beam as a continuous fluid of charged particles, 
neglecting collisions, and linearizing the one-dimensional fluid 
equations for electrons, yields 






+ V 

at + j ^ 


— E , 
m 
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Jn 3v an 

+ N . + V . = 0 , (3.2) 

3t j Bx J &x 


where N., V., n. f and v. are zeroth and first order densities and 
j j j’ j 

velocities for the j-th beam. Assuming a solution of the form 
A(x,t) = A(u),k) exp[-i(u>t - kx)] for the first order quantities, we 
have 


- i u)v . + ikV . v . s - — E , - ojn. + kv.N. + kV.n. = 0 . (3.3) 

J J J m e 3 JJ JJ 

Coupling Eq. (3.3) with the Poisson equation, 


ikE - - i-yx > <3 - 4) 

J 

yields the following dispersion relation for longitudinal oscillations 
of the beam system 
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(3.5) 




where we have assumed that the positive ions form an immobile neutral- 
ising background. 

We introduce F (V) by 

6 F j<V ' VV ' V J - V, “ « ■ f3.e> 

where the beams are spaced with equal velocity difference, 5 
Equation (3,5) can now be written as 



where is the electron plasma frequency. Dawson has shown that, 

in the limit 6 - 0 , Eq. (3.7) may be written as the sum of an integral 
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and a singular term, 



2 2 

TT U> p F(t|>/k) 

fik 2 sin 2 (ttuj/ 6k) 


(3.8) 


where denotes the Cauchy principal part of the integral, and 
F / = dF/dV . 

In order for Eq. (3.8) to have a solution, the second and third 
terms on the right hand side have to be finite when 6 — 0 . This 
requirement on the third term yields 


Im u) 



(3.9) 


which demonstrates the existence of unstable modes. Substitution of 
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this, and use of a Maxwellian velocity distribution 

F(V) = exp(-V 2 /2v^)/(2TT) 1//2 v t , in Eq . (3.8) gives finally, 

7^72 {t) 6XP (" z2 * 1 I? 1 ) = 1 + + ZZ(Z) ’ (3,10) 


(2tt) 

32 l/2 

where Z is the plasma dispersion function, with z = u>/(2 ! kv^) . 
The positive sign in Eq . (3,10) is to be used for Im u> > 0 , and the 
negative sign for Im & < 0 . 

The dispersion relation expressed by Eq, (3,10) has two complex 

conjugate roots for each value of k , Writing u >■ = a. + t 

J J J 

tr - 1A//9 1 / 2 ! 


) yields 



k6 tai" 1 

/ v” Z< V \ , 

kj6 , 

— tail 

2rr 

\l + k\ 2 + ^Re Z( ?j )/ 

= ± S{ ln 

<4n 2 v ./2 1 / 2 tt6 ) - 

X J 

(3.11) 

-f H 

'[1 + k\ 2 + Z ^ J ) ^ 2 

+ [Im Z(|J] 2 ) J • 


Each value of j denotes a particular beam, so there are clearly two 
modes for each beam. These are the normal modes of the many-beam system, 
with a Maxwellian envelope, for a given wavenumber k . Any small 
amplitude macroscopic behavior can be expressed as a weighted sum of 
these modes „ Equation (3.11) is valid only in the limit 5 0 . Note 

that B vanishes as 5 -* 0 . 

J 

Figure 3,2 is a computer simulation of the behavior of the 
streaming instability. The total electric field energy is plotted 
against time for two cases. In both, the total field energy starts 
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FIG. 3,2. Streaming instability. [L = 128 £x, h = = X D , 

6 = &v, v i = -4v t , v 2 = 4v^ , yo p At = 0.25] 

(a) § = v t /8, 8192 particles, 64 beams. 

(b) 5 = v^J 16, 16384 particles, 128 beams. 




growing almost exponentially after a certain period of time. The 
observed growth rate is about 0.09o)^ for 6 = v ^/8 * anc * about 0.05^ 
for 6 = v^/16. From these results, we may conclude that it is possible 
to carry out simulations at early times during which the field energy 
associated with wave phenomena of interest is much larger than the 
total energy of the fluctuating fields. In Fig. 3,2, for example, a 

-4 

perturbation with electrostatic energy of 10 times the thermal energy 
may be studied up to yu^t =« 40 for 6 = v ^/8 7 anc * l° n S er ^ or 
5 = v t /l6 . 

In practice, fi cannot always be made as small as is desirable, 
because the smaller 6 is to be > the more particles are necessary. 

Also, it is often necessary to follow the behavior in nonlinear problems 
for considerable periods of time. Making the beam spacing g small 
enough for such a simulation may become prohibitively expensive. As is 
shown in Section 3,4, periodic smoothing makes a long-time simulation 
possible with a relatively small number of beams. 

3,3.2 Recurrence Phenomenon 

In a simulation using the quiet start tehenique, a perturbation 
with wavenumber k first damps to a low level at the Landau damping 
rate 3 ^, and then reappears suddenly at time t = 27 x/k 6 with higher 
amplitude than its initial one. After its reappearance, the perturbation 
decays with a damping rate slightly different from the Landau damping 
rate. This is clearly demonstrated by Fig. 3.3, which was obtained from 
a result of our simulation. This recurrence results from approximating 
the continuous distribution function by a finite number of delta- 
function beams; the perturbation with wavenumber k on each beam 
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comes back into phase after a time , and with larger ampli- 

tude than its initial value due to the streaming instability. This 

g 

phenomenon was also observed by Denavit. In the limit of infinitely 
many beams, phase-mixing prevents recurrence of the initial state. 
Similar recurrence phenomena have been observed in the numerical 
solution of the Vlasov equation by the Fourier-Hermite method, 40 the 
finite difference method, 40 and the Lewis variational method. 41 



FIG. 3.3. Recurrence of initial state. Wavenumber and 
recurrence time are given by = 3n/l6 , and 

“>p T R “* 75 ‘ = V 7 ' N = 2048 ’ L = 32 AX, H = *x = \ D , 

6 = &v, v 1 = -4.5v t , v 2 = 4.5 v t> Wp &t = 0.25] . 
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The recurrence phenomenon has not been found to pose any problems 
in our simulations to be described in Sections 4-6. For example, it is 
observed that a growing perturbation does not show any irregularity at 
t = Tfl . Another example is that of a large amplitude wave, which 
involves particle trapping; as is seen in Section 4, the behavior of 
the wave demonstrates no irregularity at the recurrence time, . We 

may conclude that in strongly nonlinear cases, or in unstable situations, 
the recurrence phenomenon is not significant. It should be noted that 
periodic smoothing alone would not suffice to prevent the recurrence. 


3.4 Diffusion in Phase-Space 

Since the weighting functions mentioned in Section 3.2.2, and 
derived in the Appendix, do not conserve all of the moments of the velo- 
city distribution function when the smoothing expressed by Eq. (3.1) 
is performed, a diffusion of the distribution function occurs in phase- 

space due to reconstruction. In this section, the diffusion rate is 

6 

estimated following Denavit. 

Before smoothing, the microscopic distribution function is given 
by 


f< v > = £ V (v • V • 

j 


( 3 . 12 ) 


where f. represents the mass of particle j , with velocity v , and 
J J 

6(v) is the Dirac delta-function. The smoothed distribution function 
is 


f(v) 


•L 


E 


f w(v - V . ) , 
J J 


(3.13) 
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where w(v) is the weighting function, and A' 7 is the velocity grid 
size . 

Introducing the Fourier transforms of f(v) and f(v) 

00 

H(q) = J f(v) exp(iqv)dv = exp(iqv) , 

J 

00 

5( q ) = f f(v) exp(iqv )dv = W(q) ^ f exp(iqv.) 

J (3.14) 

where W(q) denotes the Fourier transform of the weighting function 
W (v)/Av , we obtain 

H(q) = W(q)H(q) . (3.15) 

A plot of W(q) is shown in Fig, 3.4 for the weighting functions given 
xn the Appendix. Since q represents frequency of velocity-space 
oscillation, Eq. (3.15) implies that the fine structure of the distri- 
bution function in velocity-space is smoothed out by reconstruction. 

To estimate the diffusion rate, we define 

D(q) = 1 - W(q ) . (3.16) 

Then, after m reconstructions, we have 

H m ( q ) = [1 - D(q)] m H(q) , (3.17) 

Since D(q) « 1 , for qAv/n » 1 , Eq. (3.17) may be approximated in 
this region by 

ex P[“WD (q ) ] H(q ) . (3.18) 

Thus, D (q ) represents the diffusion rate of a feature of scale 2n/q 
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W(q) 



In order to suppress the streaming instability, features of scale- 
length q =- tt/Av have to be smoothed out. Consider a perturbation 
with wavenumber k , The growth rate of this perturbation due to the 
streaming instability is obtained from Eq. (3,11) as 


P 


in IS , 
2-n A v 


(3, 19) 


when Av is small, where 6 = Av has been assumed. Assuming that the 
instability is suppressed by balancing this growth with the attenuation 
due to diffusion, we can obtain the frequency with which reconstruction 
is necessary. Thus, from 

/ k A v v \ 

“ p V- -|r in jij wlq> £ 1 ("“??)• <3 - 20) 

we have, 


t < T 
— S 


- 2jT 


k av 
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ln(l/w) 

ln(v t /Av) 


(3.21) 


where is the time interval within which at least one reconstruction 

is necessary. 

Substituting = tt/A x - tt/Xjj , v^/av = 7 * which will be used 
in most of our simulations in subsequent sections, and q ** tt/Av , 

W (q ) =-0.4 (obtained from Fig. 3.4 for the quadratic weighting function), 
into Eq. (3.20), yields 

aj t ^ 7 . (3.22) 

p s 

So far, we have only discussed diffusion in velocity-space. 
Diffusion in coordinate-sapce is treated in the same way. Replacing 
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v and q , by x and k , in Eqs . (3.16)-(3,18), indicates that there 
is attenuation of short wavelength perturbations at a rate given by 

H (k) =■ expf-DiD (k) ] H(k) , (3.23) 

111 

after m reconstructions. 



We shall now present some results of test runs with our model, 
first on the equilibrium behavior of a Maxwellian plasma, and then on 
linear wave propagation when the system is perturbed. The results will 


be compared with theoretical predictions. 


3.5.1 Equilibrium Behavior 

Our first numerical experiments were carried out on a 

Maxwellian plasma with no applied perturbation or external excitation. 

The most important parameters are the number of time- steps, , after 

which the smoothing operation is repeated, and the beam spacing, ft . 

In Fig. 3.5, the total field energy is plotted against time for 

6 = v /7 ■ ** will be seen that by increasing the frequency of 

smoothing, i.e. decreasing N , the streaming instability is suppressed; 

s 

for N g s 32 , the field energy stays roughly constant throughout the 
simulation run. To study this further, the initial energy spectrum, 
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Equilibrium behavior for various values of 
t /7, N = 2048, L = 32 £x, H = Ax = Xp > 6 = A 

4.5 V. t V =4.5 V . m At - n 9.R1 












and the time averaged energy spectrum are shown in Fig. 3.6. The 

energy spectrum for = 32 shows that mode energy tends to increase 

as time increases. For N g = 16 and 8, the mode energy seems to stay 

at roughly the same level, i.e., the streaming instability is stabilized. 

This observation agrees with the rough estimate given by Eq . (3.22). 

The total energy of the system, i.e., the particle kinetic energy plus 

the field energy, was found to be conserved to within 0.1% up to 

U> t =* 270. 

P 

In Figs. 3.7 and 3.8 are plotted analogous results to Figs. 3.5 

and 3.6 for = v /14. As expected, the growth rate of the field 

energy is greatly reduced compared with the case, 6 = v /7 . 

Figure 3.8 seems to indicate that values of N g between 16 and 32 are 

adequate to suppress the streaming instability. 

The important fact here is that for stable cases the total field 

-8 

energy is fluctuating at a level 10 times lower than the thermal 
energy of the plasma during the whole run. Since the level of fluctu- 
ations in the particle model for a system of length 32X^ with 4096 
particles is of the order of 10 *- 10 times the thermal energy, we 

have achieved 50 - 60 dB reduction in fluctuations. 

3.5.2 Linear Wave Propagation 

The purpose of this simulation was to verify predictions of 

39 

Landau damping for electron plasma waves. Waves were excited at t = 0 
by applying the perturbations 

£x. = e\ cos kx. , £v. = ev sin kx. , (3.23) 

i D i it i 

where x. , , and are the position, displacement, and velocity 
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FIG. 3.6. Time-averaged energy spectrum of the system 
shown in Fig. 3.5, 
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FIG. 3.7. Equilibrium behavior for various values of 
V [6 = V t /14, N = 4096, L = 32 AX , H = AX = X D , 
6 = &v, v 1 = -4.54 v t> v 2 = 4.54 v t , <jd At = 0,25] . 
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FIG. 3.8. Time- averaged energy spectrum of the system 
shown in Fig. 3.7, 
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perturbation of particle i , e is the amplitude, and k[= 2nn/L] is 
the wavenumber. In this simulation, only Mode 2 (n = 2, i.e., there 
are two wavelengths in the system) was excited. The results are shown 
in Fig. 3.9. For N £ 16, there is excellent agreement with the 
theoretical prediction by Langdon, 31 shown by solid lines, which takes 
into account finite-size particle effects and spatial grid effects. 

Although the fluctuation amplitude in Fig. 3.9 increases with time, 
it is important to note that the ratio of electrostatic energy to 
thermal energy, [(eE/m^^) /v“] , at t == 0 is 4.25 X 

simulation. Particle simulation with such good quality, at such a low 
electrostatic energy level, has not been feasible with previous models. 
For example, in order to reduce the fluctuation level to lo ' 6 times the 
thermal energy in a particle simulation with the same system length, 


it would require 10 3 - 10 4 times more particles than are used in this 
simulation. Since the computing cost increases roughly in proportion 
to the number of particles, it would be prohibitively expensive. In 
contrast, the computing cost in this simulation was found to be less 
than twice that with the corresponding CIC model. Suppose the smoothing 
operation is performed every 16 time-steps. One smoothing operation in 
our computer code takes about 7 sec on an IBM 360/67 for 8192 particles. 
It is equivalent to an increase of about 0.44 sec per time-step. Since 
it takes our computer code about 0.75 sec per time-step for the CIC 
model with 8192 particles and a system with 128 cells, the total com 

puting time per time-step is about 1.2 sec, 

in addition to our check on temporal Landau damping of a signal, 
we have verified the predicted linear dispersion characteristics of 
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FIG, 3,9. Simulation of Landau damping. Solid lines are 
the prediction of the Langdon theory [Eq . (2.66), £ = 0 
term only]. Dashed line is the prediction of point 
particle theory. 
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electron plasma waves. The results are shown in Fig. 3.10, and agree 

well with theory. The initial perturbations were applied for Modes 7-20 

with random phases. The electrostatic energy of the individual modes 

***■ 6 

excited was about 4 x 10 times the thermal energy at t = 0 * 

The simulation results presented in this section serve to demon- 
strate that the hybrid approach provides quantitatively accurate results 
on the collective behavior of plasma in the linear regime. Since there 
is no reason why it should not be equally effective in the nonlinear 
regime, for which analytical results are not so readily available, the 
hybrid approach is clearly a powerful tool for plasma simulation. In 
succeeding sections, we shall employ it in the study of a number of 
nonlinear problems . 
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FIG. 3.10. Linear wave dispersion. Simulation results are 


shown by circles and bars whose sizes indicate errors involved 
[6 = v t /7, N g = 16, N = 8192, L = 128 A x , H = A x = , 

6 = Av, v n = -4.5 v , v = 4.5 v , w &t = 0.251. 






4 ♦ NONLINEAR BEHAVIOR OF MONOCHROMATIC PLASMA WAVES 


4.1 Introduction 

In Section 3.5.2, we have studied numerically the collisionless 

damping to which small amplitude longitudinal electron plasma waves are 

subject. This phenomenon is due to wave-particle interaction: charged 

particles moving faster than the wave transfer energy to it, while those 

moving slower absorb wave energy, In a Maxwellian plasma, where there 

are less fast particles than slow ones in the neighborhood of the wave 

phase velocity, a net absorption of energy by the particles occurs. 

3 9 

Although first studied by Landau in 1946, the predicted damping was 

not verified in laboratory experiments until about ten years ago. 

43-46 

Spatial damping was then observed both for electron waves, and 

47 

ion waves. It was also verified that the measurements of wave disper- 
sion agree with theoretical predictions. 

Since plasma is a highly nonlinear medium, the linearized analysis 
gives only a limited description of its behavior. The question arises 
of what will happen to Landau damping and wave dispersion when the wave 

amplitude is increased . Theroetical studies of this question were first 
48 49 

made by O'Neil, and Al'tshul and Karpman. These authors found that 

the amplitude changes in time in an oscillatory manner, rather than being 

continuously attenuated. The amplitude oscillation is due to periodic 

exchange of energy between the wave and electrons trapped in the potential 

wells of the wave. The exchange occurs on a time scale of l/u) , where 

B 

1 /2 

UU [= (ekE /m ) 7 ] is the bounce frequency of an electron oscillating at 

d u e 

the bottom of a potential well, k is the wavenumber, and is the 

wave electric field amplitude. 
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The foregoing theories establish that Landau damping is valid only 
for short times, before particle trapping comes into play. The smaller 
the initial amplitude, the longer is the time for which the Landau 
solution applies. Because of analytical difficulties in dealing with 
phenomena involving particle trapping, the theoretical studies have been 
limited to special cases. For example, O’Neil assumed that the ampli- 
tude variation is so small that it may be neglected in calculating the 
48 

particle orbits. His theory is consequently valid only for 

^l/^B ^ 1 1 where y ^ is the Landau damping rate. The same restriction 

applies to the theory of Al’tshul and Karpman.^ Bailey and Denavit 

have taken into account the effects of the slowly- varying amplitude, 

and obtained essentially the same amplitude oscillation, except that 

the time at which the amplitude begins to grow again after the initial 

50 

damping is delayed. However, they still assume y « 1 , Gary 

L' B 

has treated the case > 1 analytically, and shown that the wave 

starts to decay at a rate smaller than the Landau damping rate at a 

time when the linear theory is expected to break down . 51 

The restriction on 7 ^/uig was removed in work by Sugihara and 

Kamimura, who investigated the behavior of Landau damping for a wide 

range of values . Recent work by Oei and Swanson has also 

taken into account the time-varying property of the bounce frequency . 53 

They have obtained solutions for 0 < 7/01 < 1 which are similar to 

L B 

those of Sugihara and Kamimura. An important feature of the contribu- 
tions by Sugihara and Kamimura, and Oei and Swanson r is that their theory 
is self-consistent: it includes the interaction of the electric field 

and the averaged particle velocity distribution. None of the other 
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theories mentioned so far is self-consistent. For example, O’Neil 
calculated the effect of the electric field on the particles, but not 


vice versa . 

All of the theoretical studies discussed so far assume that the 

electric field is so small that the distribution function in the 

resonant region can be expressed by a Taylor expansion about the wave 

phase velocity up to the first order term in velocity. This condition 

2 

may be written as » (v^/v^_) , where is the phase velocity 

of the wave (uu,k) . This is such a stringent condition that it is not 

easy to meet in laboratory experiments, especially when also satisfying 

the condition 7 /y) « 1 * 

L' B 

A few experimental data on Landau damping of large amplitude waves 

54 

have been reported. Malmberg and Wharton observed spatial amplitude 

48 

oscillation in qualitative agreement with the O'Neil theory modified 

55 

to fit the spatial case by Lee and Schmidt , Oei and Swanson compared 

their theoretical results with the experiments of Malmberg and Wharton, 

and found agreement on the amplitude oscillation lengths but not on the 

53 

detailed behavior of the amplitude. One of the reasons may be that 


their experimental parameters do not meet the condition 
2 

w/tDg » (v^/v^) • Specifically, they have Tg/tDg ~ 0.1 and 

2 

~ (v / v t ) for the results which exhibit amplitude oscillation, 
Franklin et al , have made detailed measurements of the spatial dependence 


of amplitude for electron plasma waves with different initial amplitudes, 
i.e,, for different values of However, for large initial ampli- 
tude, they failed to obtain results in agreement with the theory. This 


was ascribed to the appearance of sideband growth due to trapped particle 
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57 

instability. Their experimental parameters for the measured results 
corresponding to yj^ < 0.45 yield U)/^ < 4< Vp /v t ) 2 . This suggests 
that comparison of the available theories with the experiments is 
inappropriate. 

In view of the foregoing difficulties, computer simulation suggests 
itself as a means of bridging the gap between the theoretical assumptions 
and readily attainable experimental parameters. It allows conditions to 
be studied for which analytical approaches are not tractable. Such 
simulations have been carried out by Knorr, 1 using the Fourier- Fourier 
method (see Section 2.2.2), by Armstrong, 3 using the Fourier-Hermite 
method (see Section 2.2.2), and by Dawson and Shanny, 58 using the particle 
simulation model. Knorr observed a decrease in the damping rate for 
large amplitude waves at times such that Wfi t „ 1 . Armstrong considered 
the same problem and found in addition to Knorr* s results that large 
amplitude waves grow again after damping initially. He also found that 
the initial damping of a large amplitude wave is stronger than is 
predicted by the Landau theory. A similar observation of the enhanced 
initial damping was made by Dawson and Shanny. It is not certain, 
however, to what extent the large fluctuations inherent in particle 
simulation models influenced the behavior of the wave: in their com- 

putation, the field energy of the wave is of a comparable order of 
magnitude to that of the total fluctuation energy. 

One of our purposes has been to use the hybrid model described in 
Section 3 to investigate the nonlinear behavior of longitudinal mono- 
chromatic plasma waves more comprehensively than has been possible 
previously. The hybrid approach is very well suited to this study 
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because it does not generate troublesome fluctuation, and is free from 
numerical instabilities of the type encountered in the Vlasov approach 
(see Section 2,2). In Section 4.2, we shall consider amplitude oscil- 
lation and Landau damping. 

Another of our aims has been to investigate the nonlinear frequency 
shift of electron plasma waves. In a plasma of infinite extent, or of 
finite length with periodic boundary conditions, the frequency of a 
wave of large amplitude deviates from that of a small amplitude wave 
due to nonlinear effects. In an experimental plasma, in which a wave 
is excited at a fixed frequency, the shift should occur in wavelength 
instead of frequency. 

The frequency shift has been studied analytically by Manheimer 

59 60 61 62 

and Flynn, Morales and O'Neil, Dewar, and Lee and Pocobelli, 

1/2 

and found to be proportional to E^' . So far, there has been no 

report of laboratory observations of nonlinear wavelength shift for 
comparison with these theories. In Section 4.3, we shall test the 
theoretical predictions of nonlinear frequency shift against computer 
simulations carried out by use of the hybrid model, 

4.2 Amplitude Oscillation and Landau Damping 
4.2.1 Computations 

We have performed a series of computer simulations to demon- 
strate the nonlinear behavior of monochromatic electron plasma waves in 

a collisionless plasma. The electrostatic energy of the waves in these 

-4 

simulations was of the order of 10 times the thermal energy. This 
is about two orders of magnitude smaller than in the simulations of 


78 




58 

Dawson and Shanny. Some of the simulations by Knorr , 1 and Armstrong , 3 

are m our range of energy. The computations have not, however, been 

carried out for long enough times to demonstrate amplitude oscillation. 

In Fig. 4.1, we demonstrate the amplitude oscillation phenomenon 

for a large amplitude wave predicted by O'Neil . 48 m this simulation, 

8192 particles were followed in a system 64 long, divided into 64 

ceils. The continuous Maxwellian velocity distribution was replaced 

by 128 beams v t /l4 apart. Velocity-space was covered from 

-4.25 v t - 4.82 v t by a grid with mesh size equal to the beam spacing. 

Periodic smoothing was carried out after every 16 time-steps, a time- 

step being 0.25/(jjp . Mode 3 was excited initially according to 

Eq. (3.23), and the evolution of the amplitude was followed in time, 

with periodic boundary conditions applied in space, it is clear from 

Eig. 4.1 that the amplitude oscillates, as predicted. 

The initial amplitude of the wave was eE /m v m =* 3 4 v 10~ 3 

0 ' e t p * • 

corresponding to a bounce frequency of w /m ^ 0.09. The measured 

xj ir 

initial damping rate is 7^/to 0.0119. These combine to give 

V^B *“ 0 - 13 - The mea sured frequency is w / w =- 1.15. The corresponding 

X k 

wave phase velocity is Vp /v t - 3.91, so that =- 13 and 

. 2 

( V p/ v t ) =“ 15 are of the same order. 

Figure 4.2 shows the temporal behavior of the distribution function, 
in the vicinity of the wave phase velocity, averaged over space. The 
changes in the distribution are relatively small; for example, at 

V ~ 48 ’ the ratio of the peak value to the maximum value of the (nearly) 
Maxwellian distribution is of the order of 10 ~ 3 . A bump is formed in 
the distribution for u> t =- 48, and reappears for u> t =-144 and 240 

p p 
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FIG. 4.1. Amplitude oscillation of a large amplitude 

[H = Ax = X , OJ At = 0.25] . 
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monochromatic wave. 




f (v/v t ){arbitrary units) 



FIG, 4,2. Temporal behavior of the spatially averaged distribution function in the 
simulation shown in Fig, 4,1. The phase velocity of the wave is marked by an arrow 



Comparison with Fig. 4.1 indicates that these times correspond approxi- 
mately to minima in the amplitude. The height of the bump becomes 

progressi vely smaller on its reappearances , because of phase-mixing of 

48 

the trapped particles. A similar bump was observed by Armstrong, 

and considered to cause growth of waves with phase velocities lying in 

3 

that region of the bump that has positive slope. 

The bump on the tail of the distribution function has spatial 

structure. This is shown in Fig. 4.3, and may be contrasted with the 

initially spatially homogeneous distribution whose evolution is con- 

63,64 

sidered in the quasilinear theory of a warm beam-plasma system. 

The figure shows that the particles rotate by a half-cycle in phase- 

space from oOpt = 48 to 96, and another half-cycle from = 96 to 144, 

The cycle is repeated for u) t = 144 to 240. It is clear that the 

P 

phase-space structure becomes progressively less distinct as time 
increases . 

In Fig. 4.4, we present the results of a series of simulations for 
various values of the initial electric field, E^ , expressed in terms of 


the convenient parameter 7 ^/^g > where we recall that = (ekE^/m^) 

Only one mode was excited at t = 0 for each simulation run, and a 
different mode and amplitude were used in each run. The amplitude was 
normalized to unity at t = 0 in the plots. It will be seen from 
Fig. 4.4 that amplitude oscillation occurs for small values of 7 T /u) , 

and that the oscillation becomes less pronounced, with Landau damping 
extended for a longer period, as 7 /yj increases. The fluctuations 
in the curves for large values of 7 ^/^ are due to round-off 

errors made in representing numbers by a finite number of digits in 


1/2 
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FIG. 4.3. Phase-space plot at six different times, in 
the simulation shown in Fig, 4.1. The phase velocity 
of the wave is marked by an arrow. The numbers 
represent the value of the distribution function. 
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FIG. 4.4. Temporal evolution of amplitude for various initial amplitudes, i.e. 

y / 0),, , listed in Table 4.1. 
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the computer. The mode numbers and other parameters used in this 
series of simulations are tabulated in Table 4.1. 

TABLE 4.1. Parameters in the series of 
simulations presented in Fig. 4.4. 



Mode n 

L A„ 

U /u 

B / p 

7 / CO 

'll B 


a 

6 

128 

0.092 

.. . 

0. 13 

0.9 

b 

7 

128 

0.097 

... 

0.34 

1.0 

c 

7 

128 

0.061 ! 

0.55 

■B9HH 

d 

2 

32 

0.093 

0.71 


e 

2 

32 

0.079 

0.R3 

1.4 

f 

1 

15 

0 092 

0.93 

1.4 


4.2.2 Comparison with Theory 

There are a number of theories available with which we can 

31 

make comparisons; the linear theory of Langdon, the amplitude oscil- 
lation theory of O'Neil,^** and Bailey and Denavit,^ and the nonlinear 

52 

theory of Sugihara and Kamimura. 

Langdon; The theoretical values of the damping rate, y , and 
frequency, uJ T , are calculated from Eq. (2.66), which includes finite- 
size particle and spatial grid effects. Retaining only the £ = 0 
term in the summation, we obtain 7^/iV p — 0.0118 , and - 1*145 

for Mode 3 plotted in Fig. 4.1. We see very good agreement with the 
measurements described in Section 4.2.1. The theoretical predictions 
for each mode presented in Fig. 4.4 have also been found to agree with 
the measured initial damping and frequency with errors of less than 1%. 
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O’Neil: We may compare the theoretical predictions of O’Neil with 


the simulation results shown in Fig. 4.1. By solving the Vlasov equation 


for a large amplitude wave, O’Neil obtained the time-dependent damping 

, 48 

rate, 



o 2 • 

2n TT si 


Lartl 


5 2., _2n . -2n 

K F(k) (1+Q )(1+Q ) 


0 


% 2 . r<2n+l)TTtl ) 

(2n+l )tt k sxn[- F( ^ J 

+ F( K ) 2 (l+Q 2 a+ 1 )(l+Q' 2n ' 1 ) j 


(4.1) 


where F(k) [=F(k,tt/2 ) ] is the complete elliptic integral of the first 

kind, Q = exp |- n F[(l-K 2 ) 1 / 2 ]/F(K) j , and t = l/(u fi . It can be shown 

that t) , given by Eq # (4.1), reduces to the Landau damping rate in 

the limit X/t « 1 . In the time-asymptotic limit, y( t) vanishes due 

65 

to phase- mixing, and a Bernstein-Green-Kruskal (BGK) mode is formed. 

In contrast, the solution of Al'tshul and Karpman, obtained by using the 


quasilinear approximation, does not demonstrate the phase-mixing but 

49 

predicts that the amplitude continues to oscillate* O’Neil has indi- 
cated, however, that it is not certain whether their solution is correct 

. . 48 

to order y t . 

L 

We have computed Eq. (4.1), including terms up to n = 3 . We have 

substituted the numerical value y. /u) = 0.0119 obtained from our com- 

puter simulation (Fig. 4.1). For the bounce frequency, we have used 

m /(U = 0.09, calculated from the initial amplitude in the same simu- 
B' p 

lation (Fig. 4.1). The amplitude variation thus obtained is plotted 
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in Fig. 4.1. After damping initially, the wave starts to grow somewhat 

earlier than it does in the simulation. This can be ascribed to the 

change in wave amplitude, which was not taken into account by O’Neil. 

Bailey and Denavit: These authors incorporated the effects of 

. / 2 

slowly-varying wave amplitude to lowest order in a/CC , where 

a(t) = [ekE(t)/rn e ]' L / 2 , a) B = a(0) , and cc = da/dt , and obtained the 

following set of equations describing the time evolution of the ampli- 

. . 50 
tude, 


dt “ 1 - 7 ^/n ) 64n7 L U) B 


( 1 t T uT| 

2 1 R + X R J ~ 


2ttQUJ. 


B 


X cos - 2u> B t 


(S - si < 2 v>) 


u 

f 


Si(u) = / de , 


(4.2) 


T uT 

where the quantities I and I are given by 

ft ft 


i: = 


00 

E/ 


K(2a-1)(QQ 0 ) 


n-1/2 


n=1 . F(k)F(k 0 )(1+Q 211 ■Su+Qq 11 * 1 ) 


sin 


(2a-l) n f q(t / )dt " 
2 J F[K(t')] 


dK 


(4.3) 


uT 

R 


■if 


n(QQ 0 ) n 


n=l *o 


sin 




q(t / )dt 7 


"icTt' jr[K(t'yj 


dK 


J 

(4.4) 
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(4.5) 


The values of k q and n(t / ) for given values of k and t in 
Eq. (4.3) are obtained from 

a[E(K) - (1 -k 2 )F(k)] a const , 

while those in Eq. (4.4) are obtained from 

= const , (4.6) 

K 

where E(k)[ = E(k,tt/ 2)] is the complete elliptic integral of the 

second kind. We have solved Eqs . (4.2)- (4. 6) numerically, for the same 

values of y and oj used above, and with the results plotted in 

Fig, 4,1, There is very good agreement between the theory and the 

simulation. We note, however, that there is a slight difference in 

amplitude, and that the phase-mixing is somewhat slower in the simulation 

results than the theory predicts. These differences are probably due, 

. . 2 

first, to the fact that the condition, U)/U) n » (v /v ) , is not satis- 
fied in the simulation, and second, that the theory of Bailey and 
Denavit is not self-consistent. 

Sugihara and Kamimura : These authors derived from the Vlasov 

equation a set of integro-dif ferential equations which describe the 
behavior of the amplitude of a monochromatic wave. Numerical solutions 
of these integro-dif ferential equations demonstrated amplitude oscilla- 
tion for 1 > and Landau damping for 7^/^g » 1 - 

We may make a comparison between our simulation results and the 
theoretical results of Sugihara and Kamimura. First, the behavior of 
the distribution function obtained in the simulation (Fig. 4.3) may be 
compared qualitatively with that from the theory. Sugihara and Kamimura 
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presented phase-space plots for 7^/®% ~ 0*1 at three different times, 
corresponding to the first point of minimum amplitude , the first maximum, 
and to a point where the amplitude nearly ceases to oscillate (see 
Figs. 5-7 of Ref. 52). In our simulation, these times correspond to 

y$ t =* 48, 96, and 240. We find that their results and ours are consis- 

tent in the amount of rotation in phase-space, and the distinctive 
pattern. However, in our simulation, the wave amplitude still shows 
oscillatory behavior for o)^t ~ 240, in contrast to the solution of 
Sugihara and Kamimura. For > 240, our results show continuing 

particle rotation in phase-space, and a tendency to develop a circular 
plateau (see Fig. 4.3, t s 312). This is consistent with the O'Neil 
solution in the time-asymptotic limit, i.e., formation of a BGK mode. 

Next, we may compare the simulation results given in Fig. 4.4 
with those of Sugihara and Karaimura. Some of their results are reproduced 
in that figure. First, we note that their calculation shows that, for 

y /id as 0.1, the amplitude approaches a constant value after nearly two 

L B 

periods of oscillation, although the distribution function still retains 
nonuniform features. In our simulation, however, the amplitude oscilla- 
tion lasts more than two periods, and does not seem to die out so 

quickly. This fact seems to be in at least qualitative agreement with 

54 

a nonlinear spatial Landau damping experiment by Malmberg and Wharton 
in which there was no clear sign of phase-mixing , A similar feature of 
this persistent amplitude oscillation was also observed in the behavior 
of an externally excited large amplitude wave in a simulation of side- 
band instability by Denavit and Kruer. 11 Second, we recall that 
Sugihara and Karaimura found that there is a critical value of 
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y /\n = 0 77 which separates waves into those with oscillatory behavior 

r \J B 

(7 /u> < 0.77), and those which are continuously damped (Tj/iUg > 0.77). 

Figure 4.4 indicates that there is no such critical value below 
7 l/ W B = °- 93 - Third > we note that there is a tendenc y in our simulation 
results for the amplitude to decrease to a lower level, for a given 
value of 7 k <«3 > than is predicted by the theory of Sugihara and 
Kamimura; the first maximum is also lower than the theory predicts. 
Although the simulation results given in Fig. 4.4 are similar to the 
theoretical results obtained by Sugihara and Kamimura, it is important 
to note that in our simulations > ( v p / v t ) 2 ’ whereas they implicitly 

assumed that ttt/tUg » (v /v^) 


4.3 Nonlinear Frequency Shift 
4.3.1 Computations 

In Fig. 4.5, we show the variation of the nonlinear frequency 
shift of electron plasma waves as a function of the electric field 
amplitude. In this simulation, 4096 particles were followed in a system 
50 X D long, divided into 64 cells. The continuous Maxwellian velocity 
distribution was replaced by 64 beams spaced v^/7 apart. Velocity- 
space was covered from -3.79 v - 5.21 v t by a grid with mesh size 
equal to the beam spacing. Periodic smoothing was carried out every 
16 time-steps, a time-step being 0.25/tu p . Periodic boundary conditions 
were applied in space. 

Mode 3 was excited initially according to Eq. (3.23), with amplitude 
(eE g/m v t uj p ) varying from small values (9 x 10 3 ), which exhibit Landau 
damping, to large values (3 X 10 ) such as were studied in the Simula- 

c Q 

tions of Dawson and Shanny. For each simulation with different 
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FIG. 4.5, Nonlinear frequency shift of an electron plasma wave. 
[H = AX = (50/64 )\ d , o) p At = 0.25] . 



amplitude, the frequency of Mode 3 was measured by computing the total 

amount of phase change in the Fourier transform of the electric field 

between ^ t = 6 and 60, The frequency shift plotted in Fig, 4.5 was 

then obtained by subtracting the linear frequency, = 1.247 , 

obtained from Eq. (2,66), from the measured frequency. Except for very 

small amplitudes, the nonlinear frequency shift is proportional to 

1/2 

, and given by 



0.006 


0.2 — . 

(U 

P 


(4.7) 


To check the dependence of this result on the beam spacing, the simula- 
tions were repeated with the beam spacing halved, and the same number of 
smoothing operations. The differences in frequency shift were not 
more than 3%. 

A significant fact to note here is the high degree of accuracy 
with which it was possible to determine the frequency, and frequency 
shift. The model based on the hybrid approach is, therefore, much more 
efficient than a particle code in terms of computing cost for this 
measurement . 


4.3.2 Comparison with Theory 

59 

Manheimer and Flynn examined the self-consistency of the 

48 

O'Neil solution for the time-asymptotic state : they studied whether 

the potential created by the O’Neil solution satisfies the Poisson 
equation. They found that it is approximately self-consistent if a 
frequency shift given by 
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(4.8) 



is included, where p is a numerical factor equal to 2 1 / 2 , f ± s the 
tial distribution function, and is the linear plasma permittivity 

In deriving Eq. (4.8), Manheimer and Flynn only considered the trapped 
particles with simple harmonic motions, i.e., those near the potential 
wells of the wave, and the untrapped particles with straight line orbits. 
Morales and O’Neill solved an initial value problem to find the lime- 
dependent shift in the complex frequency of the wave. 60 They took into 
account the exact trajectories for both the trapped and untrapped 
particles, and obtained a frequency shift which varies in an oscillatory 
manner and approaches a constant value in the time-asymptotic limit. 

Their time-asymptotic frequency shift is expressed in the same form as 
^ ■ (4.8) except that ^ is given by 




f a J.kI2E(k) - FQpf [ 

2(e(k)-F(k)) + k 2 F(k)1 2 ) 

1 ( F(k) + 

+ g ^ > 

K F(k) 


1.63 


(4.9) 

This result is more accurate than that of Manheimer and Flynn, who 
treated particle trajectories in the approximation mentioned above. 


Lee and Pocobelli predicted frequency shifts for waves with 

V p/ V t ~ 4 up to about 50% larger than those predicted by Morales and 
O’Neil. These were obtained by including effects of electrons not in 
the vicinity of the phase velocity of the wave. 62 in contrast to these 
theories treating the case in which the wave is switched on suddenly 
at t = o , Dewar considered the case of an adiabatically excited wave 
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61 

i.e., the wave was turned on gradually. He obtained a time-asymptotic 
frequency shift similar to that expressed by Eq. (4.8), but with 
p = 1.09. 

Substituting (D^/d) = 1,247 for Mode 3, obtained from Eq. (2.66), 

and the Maxwellian distribution for f Q in Eq. (4.8), we have 


( -0.19u> b (Morales and O'Neil), 
I -0.13^ (Dewar) , 


(4.10) 


which are plotted in Fig. 4.5. We see that the slopes of the lines 
from the simulation, and from the theory of Morales and O’Neil, are very 
similar. This is to be expected because our simulation of an initial 
value problem resembles the Morales and O'Neil problem, rather than 
the Dewar problem. It should be remembered, however, that the theoreti- 
cal result is the time-asymptotic value, whereas the measured frequency 
shift is an average over the period m t = 6 to 60. It should also be 
recalled that the value of u)g corresponds to the initial amplitude 
of the wave. Since the theoretical result due to Morales and O'Neil 
was obtained under the condition that the amplitude variation is very 
small, it does not matter much whether the bounce frequency is computed 
from the initial amplitude or from the time-asymptotic amplitude. In 
our simulation, however, the amplitude variation is not negligible; if 
the bounce frequency were computed from the time- asymptotic amplitudes, 
the points in Fig. 4.5 would be moved towards the theoretical line of 
Morales and O’Neil. 
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4.4 Summary 


The nonlinear behavior of monochromatic plasma waves has been 
studied over a wide range of wave amplitudes, by use of the low-noise 
model based on the hybrid approach. 

In the study of amplitude oscillation and Landau damping in 
Section 4.2, we have attempted investigation in areas where analytical 
approaches are not easily tractable, i.e,, in cases where the condition, 
» (v /v^) , is not. satisfied. The results of our simulations show 
good qualitative agreement with the theories of Bailey and Denavit, 



However, there are significant differences between our simulation results 

and the theoretical results of these authors; first, phase-mixing of the 

amplitude oscillation is slower than predicted, and second, there exists 

no critical value of y An within our parameter range such as was 

w B 

found by Sugihara and Kamimura. These results will be helpful in better 

understanding the phenomenon, and in developing an analytical theory in 

2 

cases where <* (v^/v^) 

In the study of nonlinear frequency shift, in Section 4.3, we have 
measured the frequency shift for finite amplitude waves, and compared the 
results with theoretical predictions. It has been demonstrated that the 
simulation results agree well with the theoretical predictions of Morales 
and O’Neil. 
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5. SIDEBAND INSTABILITY 


5.1 Introduction 


In 1968, Wharton et al, reported results of experiments on large 
amplitude, longitudinal electron plasma waves, in which not only was 
spatial amplitude oscillation of the type discussed in Section 4 
observed, but also spatial growth of sidebands separated from the 
frequency of the large amplitude wave by the bounce frequency of trapped 


electrons , 


These experiments have stimulated a number of 


theoretical studies that may be classified into two types of approaches. 

One is based on a wave-wave interaction mechanism between the large 

57 67 

amplitude wave and sideband waves (Kruer et al.; Goldman; Goldman 

68 69 70 71 

and Berk; Wong; and Mima and Nishikawa ' ). The other is a quasi- 

linear approach based on wave-particle interaction (Shapiro and 

72 73 74 75 

Shevchenko; Bud'ko et al.; Manheimer; Yagishita and Ichikawa; 

7 6 

and Brinca ). 

Other laboratory experiments on sideband instability have been 

77 78 

carried out by Franklin et al., and Jahns and Van Hoven for electron 

79 

plasma waves, and by Ikezi et jal. for ion waves. These experimental 

results have verified some of the predictions of wave-wave interaction 

77 79 

theory applied to a spatial case. 9 However, there are some obser- 
vations which suggest the quasilinear mechanism as an alternative cause 

78 

of the instability. 

Computer simulations of the sideband instability have been performed 

80 11 8l 

by Kruer and Dawson, Denavit and Kruer, and Rosen et al. Kruer 

and Dawson studied the instability in a one-dimensional plasma driven 


by an external electric field of a given frequency by use of a particle 
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simulation model. They observed the growth of sidebands having a 
frequency separation consistent with the experimental results of 
Wharton et al, It was demonstrated later by Rosen et al . that artificial 
removal of trapped particles eliminates the sideband growth. Denavit 
and Kruer carried out simulations of a similar problem to make compari- 
son between the particle simulation and the Vlasov approach, and found 
close agreement between the results of the two approaches. In these 
simulations, electrostatic energy of the large amplitude wave was 
0. 1-1.0 times the initial thermal energy of the plasma. The high wave 
energy is required in the particle code since such simulations are much 
noisier than real plasmas, as discussed in Section 2. 

In laboratory experiments, however, electrostatic energy of the 
large amplitude wave is typically much smaller (10~ 4 - 10 _3 ) than the 
thermal energy. Also, in most of the analytical approaches mentioned 
above, it is assumed that the wave energy is much smaller than the 
thermal energy. In what follows we shall investigate the sideband 
instability, using the simulation model based on the hybrid approach, 
in the parameter range which allows us both to make more quantitative 
comparisons with existing theories, and to model more satisfactorily 
the conditions appropriate to laboratory experiments. The quantitatively 
accurate results obtainable by use of this model should lead to better 
understanding of the phenomena involved, and refinement of the theory. 

5.2 Computations 

In the series of simulations to be described in this section, we 
have considered an initial value problem, and imposed perturbations at 
time t = 0 according to Eq . (3.23). Mode 13 was chosen for the large 
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amplitude wave, which we will refer to from now on as the ’main* wave. 

Five simulation runs (A - E) were made with the initial main wave 

amplitude in the range 0,06 < eE^/m^v^u^ <0.5, In terms of elect ro- 

-3 -i 

static energy, this is between 1.8 x 10 and 1-2 x 10 times the 
thermal energy. Waves were also excited initially as sidebands of 
Mode 13 according to Eq. (3.23), but with random phases at the energy 
level of 10 times the thermal energy. 

In all of the computations, 16384 particles were followed in a 
system 256 long, divided into 256 cells. The continuous Maxwellian 
velocity distribution was replaced by 64 beams spaced v^/7 apart. 
Veloci ty-space was covered from -3,79 - 5,21 v^ by a grid with mesh 

size equal to the beam spacing. Periodic smoothing was carried out 
every 16 time-steps, a time-step being 0.25/uo^ . Periodic boundary 
conditions were applied in space. 

Results of a typical simulation are given in Fig. 5.1. It shows 
the evolution of the main wave and two test waves. The initial ampli- 
tude of the main wave was eE /m v j 0.12, corresponding to 

0' e t p 

l/2 

U) /& =“0.19 y where m [= (ek E /m ) ' ] is the bounce frequency of an 

Jj P d U U G 

electron at the bottom of a potential well of the main wave with wave- 
number . Mode 11 decays first, and then begins to grow. Mode 12 

shows an evolutionary pattern similar to that of Mode 11 at early times, 
and then begins to grow more slowly than Mode 11. The temporal behavior 
of the velocity distribution function at the early times is shown in 
Fig. 5.2. It is a plot of the spatially averaged distribution func- 
tion in the vicinity of the phase velocity of the main wave. The 
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RELATIVE MODE ENERGY 




FIG. 5.1. Temporal evolution of sideband instability: 

(Initial main wave electrostatic energy/ thermal energy) = 7.42 x 10 




f (v/vj ) (arbitrary units) 



FIG. 5,2. Temporal behavior of the spatially averaged 
distribution function at short times in the simulation 
shown in Fig. 5.1, Phase velocities of Modes 11-13 are 
shown by arrows, w* is the velocity given by Eq. (5,13). 
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distribution function at later times, when the sidebands are growing 

steadily, is shown in Fig. 5.3 as a phase-space plot. It is seen 

clearly that there are particles trapped near the bottom of the wave 

potential wells. We have confirmed that phase-space structure similar 

to that shown in Fig, 5.3 persisted throughout the simulation run. This 

may be compared with Figs. 3 and 8 of Ref. 80, which show only a 

negligible number of such particles to be present. In their simulation, 

using a particle code, the number of particles with the velocities in 

the vicinity of the phase velocity 4 v ) of the main wave is very 

small, because the Maxwellian velocity distribution assumed falls off 
2 2 

as exp (-v /2v t ). Since it is these particles that will stay trapped 
near the bottom of the potential wells, a situation such as observed 
in Ref. 80 may occur. 

In Fig. 5.4 are given the results of two simulations with a larger 
amplitude of the main wave than the previous simulation. The growth of 
Mode 11, in the lower sideband (k < k Q ) , is shown in Fig. 5.4(a) both 
for growth from noise, and when it is excited as a test wave at t = 0 , 
Figure 5.4(b) shows growth of a test wave in Mode 15, in the upper side- 
band (k > k^). Though not shown, we also observed growth from noise 
of Mode 15, and other modes, in the simulation carried out without test 
waves. The distribution function demonstrated similar behavior to that 
shown in Figs. 5.2 and 5.3. 

An energy spectrum obtained in one of the simulations is shown in 

Fig. 5.5. It was measured at time <*> t = 175, when the main wave 

P 

amplitude reaches nearly the fifth minimum in the temporal evolution. 

The lower sideband is seen to be larger than the upper sideband. It 
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FIG . 5.3. Phase-space plot at uj^t = 96, in the 
simulation shown in Fig. 5.1. The phase 
velocity of the main wave is marked by an arrow. 
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RELATIVE MODE ENERGY 



FIG. 5.4(a). Temporal evolution of sideband instability: 

(Initial main wave electrostatic energy/thermal 

^2 

energy) = 2.97 x 10 . Main and lower sideband waves. 
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RELATIVE MODE ENERGY 



FIG. 5.4(b). Temporal evolution of sideband instability: 

(Initial main wave electrostatic energy/ thermal 

-2 

energy) = 2.97 x 10 . Main and upper sideband waves. 
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RELATIVE MODE ENERGY 



FIG. 5.5, Energy spectrum at u^t = 175 due to sideband 

instability: (Initial main wave electrostatic energy/ 

-2 

thermal energy) = 2,97 x 10 
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should be remembered, however, that the test waves initially decay and 
then begin to grow. Consequently, the fastest growing mode (Mode 10) 
does not necessarily appear as the sideband peak. Some of the modes, 
for example Mode 14, are still below their initial level at the time 
when the energy spectrum is measured. 

A case with a heavily-damped main wave is presented in Fig. 5.6. 
in this simulation, Mode 17 was used as the main wave. The test wave 
shown was the fastest growing mode. We see that the main wave first 
decays, and then undergoes amplitude oscillation with slow damping. The 
test wave shows continued growth after it has reached an energy level 
comparable to, or exceeding, that of the main wave. 

From the results shown, we note the following: 

(a) The main wave amplitude shows oscillatory behavior, 
corresponding to that predicted by O’Neil for the temporal case. 

(b) Test waves exhibit initial decay which is stronger than the 
corresponding linear Landau damping, and then show approximately 
exponential growth, 

(c) Both sidebands grow from noise with growth rates corresponding 
to those of the test waves. 

(d) The lower sideband is higher in amplitude than the upper 
sideband . 

(e) There is some modulation superimposed on the growth of the 
sidebands. This seems to be correlated with the amplitude oscilla- 
tion of the main wave. 

(f) In the case of a heavily-damped main wave, test waves exhibit 
features similar to (a) through (e). 
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RELATIVE MODE ENERGY 



PIG. 5.6. Sideband instability with a heavily-damped main wave: (Initial 

main wave electrostatic energy/thermal energy) = 7.1 * 10~ 3 . The damping 

rate and phase velocity of the main wave are y /yj =* o 085 and v /v =-3 1 

L p ‘ ' p / t 

[Eq. (2.66)]. 



The foregoing features are common to all of our simulations, and 
suggest that it might be profitable to divide theoretical description 
into two parts; first, the transient processes occurring at the 
earliest stages of evolution, and second, the development at later times 
when the sidebands are growing steadily. 

5 , 3 Comparison with Theory 

^ 57 

In order to explain the sideband instability, Kruer et al , 
considered a simple theoretical model in which the trapped electrons are 
treated as a bunched beam of harmonic oscillators of frequency , 

and obtained a sideband growth rate consistent with the experimental 

67 

results of Wharton et al. Goldman took a rather different approach. 

.48 

It is known that the time-asymptotic limit of the O’Neil solution 

(2 C 

is a large amplitude BGK mode. As long ago as 1962, Pfirsch had 
speculated that a large class of BGK modes might be unstable, although 

QO 

he did not pursue the question. Goldman examined the stability of the 

BGK modes, and showed that the sideband growth is due to a parametric 

type of coupling between waves enhanced by the trapped particles. He 

obtained the results of Kruer et al. as a special case in which trapped 

electrons are localized at the bottom of the potential wells of a BGK 

wave. Goldman and Berk obtained a dispersion relation, in the bunched 

beam approximation, including the contribution of trapped electrons to 

68 

the frequency shift of the large amplitude wave. They showed that 

69 

the growth rate is enhanced above that of Kruer e t al . Wong has 
investigated the stability of two types of BGK modes, including the 
effects of resonant interaction with both trapped and untrapped 
electrons, and obtained results similar to those of Kruer et al . Mima 
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and Nishikawa have developed a theory without assuming any particular 

form of BGK mode, and predicted sideband instabilities over two wave- 

number ranges given by |k - k Q | = (2n + l) 1 / 2 W B /v p or 

jk - k Q j « a) B /v p , where v^ is the phase velocity of the large 

70 

amplitude wave, and n(^0) is an integer. Mima and Nishikawa 

later investigated the stability of a BGK mode whose untrapped particle 

distribution is chosen to be that given by the Landau linear theory. 71 

All of the theories discussed so far concentrated on investigation 

of the wave-wave interaction mechanism. An alternative approach is a 

79 

quasilinear wave-particle interaction. Shapiro and Shevchenko, ^ and 
73 

Bud’ko et al., have studied the excitation of sidebands due to 

resonant wave-particle interaction, using the O’Neil time-asymptotic 
48 

solution as their starting point. Bud’ko et al . found that only a 

lower sideband satisfying [u - kv I =* 0.9y) can be excited. However, 

for the parameters used in the experiment of Wharton et al . , their theory 

does not predict instability. Shapiro and Shevchenko used a different 

distribution function for untrapped particles, and found that both 

sidebands can be excited with different growth rates. Substitution of 

the experimental parameters of Wharton et al . in their theory yields 

growing solutions, but the growth rates are much smaller than those 

57 

obtained by Kruer et al. The stability of the O’Neil solution has 

also been investigated by Manheimer, who considered only the particles 

trapped near the bottom of the potential wells, and predicted that the 

74 

lower sideband is unstable while the upper sideband is stable. His 
theory is considered to be a lower order approximation to the theories 
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of Shapiro and Shevchenko, and Bud’ko et al. In the quasilinear theory 

75 

of Yagishita and Ichikawa, the stability of the time-asymptotic distri- 

49 

bution function found by Al*tshul and Karpman (see Section 4) was 

studied. It was shown that the trapped electrons can cause sideband 

instability through interaction with externally excited test waves. 

The theories mentioned so far examine the stability of either a 

stationary or time- asymptotic equilibrium state, involving a large 

amplitude wave. In contrast, Brinca has used the quasilinear theory 

during the transient following application of a large amplitude wave 
76 

at t = 0. He determined the variation of the sideband growth rate, 

as a function of time, from the slope of the averaged velocity distri- 
bution function in the vicinity of the phase velocity of the large 
amplitude wave. 

Examining the theories discussed so far, in the light of the results 
of our computer simulations in Section 5.2, it seems appropriate to 
compare the simulation results at early stages of temporal evolution 
with the theory of Brinca, and at later times with the theory of Kruer 
et al . In what follows, we shall, first, consider the quasilinear 
wave-particle interaction theory due to Brinca, and then the wave-wave 
interaction theory in the bunched beam approximation due to Kruer et al. 

5.3,1 Quas i 1 i near Theory 

Theory : Consider a large amplitude, electron plasma wave, 

excited at time t = 0 , in an infinite collisionless plasma. The 
electron motions are described by 
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dx dv e 

dt = v ’ at = " IT E o sin <V ’ k o x) ' (5 - 1} 

e 

where x and v are the position and velocity of an electron, and 
is the frequency of the large amplitude wave. Since E^ is 
assumed to be constant, Eq. (5.1) yields the first integral of the 
motion, 


W 




eE 

+ cos (m 0 t - k Q x) . 

0 


(5.2) 


The distribution function in the presence of the large amplitude 
wave is determined by the Vlasov equation, 


9£ + v-E 

at + ^ 



sin (u) 0 t - k Q x) 



(5.3) 


Since Eq . (5.1) represents the characteristics of Eq . (5.3) in phase- 
space, the solution of Eq. (5,3) is given by 

toc.v.t) => f 0 + 2- ; 0 (*.».t) t' , (5.4) 

in the resonant region, v =* . Here, is the initial electron 

velocity distribution function, f^ denotes the derivative of 

with respect to v , and 5 0 (x,v,t) is given by solution of Eq. (5,1) as 


5o - *7f dn [ P<K '5 ) T l;] ' <« 2 <1). 

tor unt rapped particles, and by 

cn [ F (\ > 0* (fi2 > 1} - 


(5.5) 


(5.6) 
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for trapped particles, where F represents the elliptic integral of 
the first kind, dn and cn are Jacobian elliptic functions, 
denotes the initial value of g , and k , t , § and £ are defined 
by 

» \V 2 


2 2eE 0 
K = 


k o w+eE o 




(5.7) 


1 = 2 < k 0 x - - tt) , sin £ = k sin g . 

The positive or negative sign is used in Eqs. (5.5) and (5,6) according 
to whether £ , § > 0 , or £ , g < 0 . Equation (5.4) is valid pro- 
vided that the amplitude is small enough to allow a Taylor expansion 

of the distribution function about the phase velocity of the large 

/ 2 

amplitude wave, i.e. ^ B /^ 0 « ( k Q V t' x o^ ‘ 

Having found the solution of the Vlasov equation, the next step 
is to average Eq. (5.4) so as to remove fast oscillations occurring 
on a time scale l/uo^ anc * on a l ea gth scale l/k^ > and to obtain a 
slowly-varying velocity distribution (f(v,yj t)) . After some mani- 

pulation, the averaged solution 




(5.8) 


is obtained, in which 


TWI,, 

£w > = 2 

0 k Q KF(K) 



-.2n 


2n 2 
(1+Q ) 


cos 


r n ™vil 
L J ) ’ 


(5.9) 


(w) = w* 


E(k) 


(K < 1) , 


(5.10) 
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for the untrapped particles, and 


<* ft > « 


8tTU) 


B 


O' k 0 F(T]) 


n=l 


Q 2 * 1 " 1 

(1 + Q 2n_1 ) 2 


cos 


r<2n-l)m B t-l I 
L 2F(T|) J ( ’ 


(5.11) 


<w) = w* [E(Tp - (1 - Tj 2 )F(T|) ] (t| 2 = ^ < l) , 

\ k ' 


(5.12) 


for the trapped particles. In Eqs . (5 . 9)-(5 . 12 ) , F (*)[:= F(k,tt/ 2)] 
and E(k)[= E(k ? tt/ 2)] are the complete elliptic integrals of the first 
and second kinds, and 


w 


Q 



W* 2= 


4uj. 

Trk, 


B 


exp | -nF[ (1 -k 2 ) 1//2 ]/F(k) j . 


(5,13) 


Since k: and 7] are related to the mean velocities defined by 

Eqs. (5.10) and (5,12), we see that ( W Q ) > which is a function of k 
or , is itself a function of the mean velocity. 

In Fig. 5.7, ( w q) is plotted as a function of velocity, with time 
as a parameter. It will be recognized from Eq . (5.8) that a plot of 
(w 0 ) indicates the shape of the averaged velocity distribution function 
in the resonant region. It will be noted that ( w q) develops finer 
and finer structure as time increases, resulting in progressively larger 
local slope of the velocity distribution function. The application of 
the theory is limited to the transient process at the initial stage of 

evolution, before the fine structure develops, i.e. up to u)t ^ 2rr . 

B 

If the averaged distribution function given by Eq. (5.8) is 
considered as a slowly-varying * equilibrium r distribution function, 
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FIG. 5.7, Temporal behavior of the averaged distribution function in the resonant region, 
the phase velocity of the main wave. T1 
test waves, (Adapted from Fig. 1 of Ref. 76.). 


v q is the phase velocity of the main wave. The other v^ f s are phase velocities of 



then a small perturbation, , due to application of a test wave may 

be described by the linearized Vlasov equation in the form 


Sf l af l e 

— i + v — - - — E = 0 

at ^ m e i av 


a<f> 


(5.14) 


where E^ is the perturbed electric field. Applying linear stability 
analysis for perturbations of the form exp[-i(yjt-kx>] , propagating 
with phase velocities near that of the large amplitude wave, yields 
a growth rate of 


>■<*> =l(r)(“ -s)(^) w/k - (5 - 15) 

where y « u) has been assumed. 

Having obtained this result, Brinca simplified it by assuming that 
the test waves are described by the warm plasma dispersion relation 


e p (k,u>) 


1 - 


U) 

p 

2 2 2 
U) -3k v t 


= 0 


(5.16) 


Equations (5.8), (5.15), and (5.16) are then sufficient to obtain the 
time evolution of the growth (or decay) rates of the test waves. 

Figure 5.8 presents some calculations. It will be seen that the side- 
bands decay initially (if they are above the thermal fluctuation level), 
and then start to grow; the larger the phase velocity separation between 
the test wave and the large amplitude wave, the stronger the initial 
damping. 
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Comparison with C omputations: The theoretical growth rate for the 


test waves in the initial development stage has been computed using 

Eqs . (5.8), (5.11), (5.12), (5.15) and (5.16). Equations (5.11) and 

(5.12) are used since |u)/k - oj^/k^| < w* holds for Modes 11 and 12 as 

shown in Fig. 5.2. The calculated results are shown in Fig. 5.9 for 

comparison with the growth rates measured from the simulation results 

presented in Fig. 5.1. Although the evolutionary patterns for Mode 12 

from the theory and computation 2 'esemble each other to some extent, the 

agreement is not good. It should be borne in mind that the theory is 

applicable only to cases in which the main wave amplitude variation is 

negligible; in our simulation the amplitude actually varies by a 

factor of more than two during the period up to u> t = 30 . To obtain 

P 

better agreement, we have used information from the detailed plot of 

the averaged velocity distribution function in Fig, 5,2. We have 

computed the growth rate, using Eq . (5.15), at intervals of cu t = 4 

P 

from the local slope of the averaged velocity distribution function at 

the phase velocities of the test waves. The results are shown in 

Fig, 5.9. There is a striking similarity in evolutionary pattern for 

both modes between the measured growth rate and the theoretical one; 

there is a difference in a; t of about 3 for both modes until 

P 

10 , and a difference of about 15 for Mode 11, and 10 for Mode 12, 
thereafter. 

In seeking an explanation for this phenomenon, it should be recalled 
that in the quasilinear theory of Brinca there is an implicit assumption 
that the waves respond instantaneously to the slope of the averaged 
velocity distribution function. However, since transient phenomena 
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FIG. 5.9. Test wave growth rate at short times in the simulation shown in Fig. 5.1. 
Brinca (1) is the theoretical result calculated by the use of warm plasma dispersion 
relation [Eq. (5.16)]. Brinca (2) is the theoretical result obtained from Fig. 5.2. 







are involved, it is more reasonable to assume that some time elapses 
before the macroscopic effects of the resonant wave- particle inter- 
action appear. We may estimate this time delay as follows. From 
Fig. 5.2, the width of the bump in the velocity distribution function 
may be estimated to be fiv =* 0,6 . Considering this bump as a set 

of streams with continuously distributed velocities, providing a con- 
tinuous range of frequencies kv , we obtain the rise time, , of 

a perturbation from the approximate equality 

k&VT . (5.17) 

d 2 

Substituting appropriate numerical values into this expression yields 

U> t , ^ 10 for Mode 11, and 9 for Mode 12, These results do not fully 
p d 

account for the discrepancies, of course, but provide a good intuitive 
explanation. 

5.3,2 Wave-Wave Interaction Theory 

Theory: The simplest model of sideband instability which 

incorporates wave-wave interaction is that involving the bunched beam 
approximation, In this approximation, the equilibrium distribution 
function is assumed to contain trapped particles localized at the 
bottom of the potential wells of a large amplitude wave propagating at 
phase velocity v^ . The trapped particles act coherently as harmonic 
oscillators of frequency Ui fi [= (ek Eg/mp 1 / 2 ] . The effects of trapped 
particles other than those near the bottom are neglected. 
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The trapped electrons are governed by the equation of motion 


—T = - %<V x n0“ V p t) ' ~ L 2 / r E(k,tu)exp[- i (wt-kx^ ]dkdu) , 
dt (2rr) J e 

(5.18) 


where x is the position of a particle in the n-th potential well, 
n 

x ^ + v t is the position of the n-th well, and E(k,u)) is the Fourier 
nO p 

transform of the perturbation electric field . The effects of the 
electric field of the large amplitude wave are taken care of by the 
first term on the right-hand side of Eq * (5.18). 

Treating the trapped electrons as a source charge density intro- 
duced into a plasma of permittivity e (k,u)) , we have 

P 

ike (k, uj)E(k,uj) = P^’^ . , (5.19) 

P € 0 


where p(k,uO) is obtained from the Fourier transform of the displace- 
ment of the trapped particles given by Eq . (5.18). Some manipulation 
of Eqs . (5.18) and (5.19) yields 


(JO, 


E (k, oj) = 


T 


. u ,2 2 
(m-kv ) -|H 

" p/ wg n 


E(k+nk Q) (j>fnu) o ) 

e (k,«j) 

P 


v 2 V /2 


(5.20) 


W T ~ 1 e m x 

0 e Oi 


x - in 

*0 “ k 

0 


- k v , 
Op' 


whei'e is the number of trapped electrons in each potential well. 

Equation (5.20) shows that perturbations at w,k are coupled to an 
infinity of perturbations at aH-n^Q , k+nk^ . 
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Since plasma does not support wave propagation at frequencies 
greatly different from the plasma frequency, , the two waves E(k t tu) 

and E(k-2k^, (U“2^) may be expected to be dominant for cu ~ (ju^ * 
Retaining only these two terms yields two coupled mode equations for 
E(k,tw) and E(k~2k^ J u)-2u>q) . The dispersion relation results from 
equating the determinant of their coefficients to zero. We obtain 



where Q = ^ - kv . If the large amplitude wave is not too large, the 
P 

warm plasma approximation for e^Ck,^) [Eq. (5.16)] may be used. 

Comparison with Computations : In Fig. 5.10 are plotted the theo- 

retical growth rates for Modes 8-12, obtained by solving Eq , (5.21) 
combined with Eq. (5.16). To make quantitative comparisons, the growth 
rates of the unstable modes were measured in the five simulations (A-E) 
described in Section 5.2. The growth rates were obtained from energy/ 
time plots for each unstable mode, similar to those shown in Figs. 5.1 
and 5.4, and plotted in Fig. 5.10. The errors involved in measuring 
the growth rates of the sidebands were about 10%. We see that there is 

good agreement for small values of yj , but that the measured growth 

B 

rates tend to be larger than the theory predicts for large values of 


UUg . In considering the discrepancies, it should be remembered that 

the theory is not valid for very large amplitude main waves, i.e. we 

require E^/4nnk T « 1 , and that values of the bounce frequency, m , 
0 7 B B 

used in this plot are not those corresponding to the initial amplitudes 
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FIG. 5.10. Sideband growth rate at later times vs. bounce 
frequency. Solid lines are solutions of Eq. (5.21) combined 
with Eq. (5.16). 



of the main wave, but have been estimated from the expected time- 
asymptotic amplitudes indicated by dashed lines in Figs. 5*1 and 5,4. 

5 . 4 Comparison with Experiments 

In comparing the results of our simulations with those of labora- 
66 77 78 

tory experiments, * ' it should be noted that the simulations were 

carried out for an initial value problem, rather than a boundary value 

problem. Since laboratory experiments deal with spatial phenomena, 

direct quantitative comparison may not be appropriate. Under certain 

conditions, however, it may be possible to transform an initial value 

problem into a boundary value problem so as to allow quantitative 

comparison with the experimental data* For example, as found by Lee 
55 

and Schmidt, the O’Neil solution for a temporal case can be trans- 
formed into a spatial solution by replacing the normalized time, 
u^t , and the parameter 7 h /% , by x/\ and , where = w/km 

and p L = 7 l /(& . 

In making this transformation, it should be remembered that in the 
simulation the system length is finite, and only a finite number of 
wavenumbers are available with equal separation 2tt/L . As a conse- 
quence, it may well be that the fastest growing mode observed in the 
simulation is not the fastest growing mode predicted by theory for an 
infinite plasma. This implies that an accurate measurement of the 
dependence of the sideband peak frequency separation, {#& , and side- 
band growth rate, y , is not available from our simulations. Neverthe- 
less, the data plotted in Fig. 5.10 seem to suggest y cc , with 

a > 1/2 , which is to be compared with the experimental observations 
that <x , 7 oc E^/^ by Franklin et al *, 77 that 
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l/2 B l.a (7 is a constant) by Jahns and Van Hoven, 

* * V ' ? ' 7 ® “ ° V 66 These measurements could be 

hv Wharton. et al* 
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amplitude is relatively small (ujg/u)p < 0,2), and the test wave growth 
(or damping) rate is not too large compared with its frequency 
(7 « tu) . 

The later stages of evolution of the sideband instability, after 
one or two phase-space rotations of the particles have been completed, 
have been shown to be well described by the wave-wave interaction theory 
in the bunched beam approximation of Kruer et al. For very large main 
wave amplitude > 0.2), the simulation gives higher growth rate 

than the theory predicts. 

Comparison with the laboratory experiments has shown that many 
features of the instability observed in our simulation at similar signal 
levels are consistent with the experimental observations, account always 
being taken of the fact that the simulation is for temporal evolution, 
and laboratory experiments for spatial evolution. The dependence of 
the frequency separation, , and sideband growth rate, y , upon the 
wave amplitude, E Q , still need to be checked. 
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6. SATELLITE GROWTH 


6,1 Introduction 

The computer simulations to be described in this Section were 

stimulated by laboratory observations by Jahns and Van Hoven of 

satellite growth at frequency & 9 (= ) when a large amplitude 

signal at , and a small amplitude signal at were excited 

83 

simultaneously. Jahns and Van Hoven interpreted the satellite growth 

as being due to four-wave passive coupling. They applied a perturbation 
84 

expansion method of solving the Vlasov and Poisson equations to a 

spatial problem appropriate to their experiment, and obtained a solution, 

85 

describing the spatial evolution of the satellite. However, the 
predicted dependence of the satellite growth rate upon the amplitude of 
the signal at did not fit the observed dependence. Jahns and 

Van Hoven ascribed the discrepancies to dissipation, higher-order 
processes, and wavenumber mismatch. 

DeNeef made similar observations to those of Jahns and Van Hoven 
in his experiments with a large amplitude wave and a small amplitude 

Q G 

wave launched simultaneously. He considered the small amplitude wave 
as a slow modulation of the amplitude and phase of the large amplitude 
wave, and calculated the amplitudes of the small amplitude wave and 
the satellite wave as a function of position. His calculation showed 
agreement with the experiments for the former, but not for the latter. 

In particular, the energy level of the satellite wave observed in the 
experiment was 90 dB below the theoretical prediction, DeNeef 
suggested that the discrepancy might be due to the strong dependence 
of the satellite behavior on the nonlinear wavelength shift of the 
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large amplitude wave. In his theory, he used the wavelength shift 

60 

calculated from the theory of Morales and O’Neil, 

Brinca considered such a process for the analogous temporal problem 
i 11 which the synchronism relations 

2k () = k x + k 2 , 2u + m 2 , (6.1) 

87 

are satisfied. He obtained coupled-mode equations which describe 
the temporal evolution of the wave amplitudes. The theory failed, 
however, to give either the observed rapid growth rate, or the observed 
satellite energy level . 

In what follows, we shall demonstrate good agreement between com- 
puter simulations and theoretical predictions based on DeNeef’s method 
applied to a temporal problem rather than a spatial one. Unlike DeNeef, 
we use the measured nonlinear frequency shift in the calculation of the 
wave evolution. The computer simulation is discussed in Section 6.2. 

The theory is described in Section 6.3, and compared with the simulation 
in Section 6.4. 

6.2 Computations 

In the simulations to be described, 16348 particles were followed 

in a system 256 \ D long, divided into 256 cells. The continuous 

Maxwellian distribution was replaced by 64 beams spaced v^/7 apart. 

Velocity-space was covered from -4.5 v^ - 4.5 v^ by a grid with mesh 

size equal to the beam spacing. Periodic smoothing was carried out 

every 16 time-steps, a time-step being 0.25/uu . Periodic boundary 

P 

conditions were applied in space. 
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Perturbations of the form given by Eq. (3.23) were applied in 
Modes 13 and 12 at t = O' . Figure 6.1 shows the results of the simu- 
lation for two different main (large amplitude) wave amplitudes, and 
three different test (small amplitude) wave amplitudes. It will be seen 
that in each case the main wave evolves almost exactly as a single 
large amplitude wave (compare with Fig. 4.1), i.e. the main wave ampli- 
tude is not large enough to cause appreciable sideband growth, due to 
trapped particle instability of the type studied in Section 5, on the 
time scale for which the satellite grows and saturates. The test wave 
follows a very similar evolution to that of the main wave. The satellite 
grows from noise, saturates at a level comparable to that of the test 
wave, and seems to show oscillatory behavior thereafter. 

From Fig. 6.1(a)-(c), we observe that the behavior of the test 
wave and the satellite is almost identical for different test wave 
amplitudes, except that the curves are shifted vertically by an 
amount which scales linearly with the test wave amplitude. This is so 
only when the test wave amplitude remains small compared with the main 
wave amplitude. If increased progressively it finally disrupts the 
particle trapping by the main wave, and hence affects the main and test 
wave evolution. From Fig. 6. 1(c)- (e) it will be seen that the growth 

rate of the satellite seems to decrease as the main wave amplitude 

. 88 

deci'eases . 

When the roles of Modes 12 and 14 are switched, Mode 12 is observed 
to grow from noise, reach the level of Mode 14, and finally exceed it. 

We have confirmed that halving the beam spacing changes the results only 
in minor details . 
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RELATIVE MODE ENERGY 



FIG. 6.1. Temporal evolution of main, test, and satellite waves, 
(a) (Main wave electrostatic energy/thermal energy) = 1.86 x 10 
(Test wave electrostatic energy/thermal energy) = 4.18 x 10 
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FIG. 6,1* Temporal evolution of main, test, and satellite waves , 
(b) (Main wave electrostatic energy/thermal energy) s 1,86 x 10 
(Test wave electrostatic energy/thermal energy) = 1,16 x 10 


L3 



RELATIVE MODE ENERGY 



FIG. 6*1, Temporal evolution of main, test, and satellite waves, 
(c) (Main wave electrostatic energy/thermal energy) = 1,86 x 10 
(Test wave electrostatic energy/ thermal energy) = 1,16 x 10 


131 



RELATIVE MODE ENERGY 




FIG. 6.1. Temporal evolution of main, test, and satellite waves, 
(d) (Main wave electrostatic energy/ thermal energy) = 7.25 x 10 
(Test wave electrostatic energy/ thermal energy) = 1.16 x 10 



RELATIVE MODE ENERGY 


I 



FIG, 6,1. Temporal evolution of main, test, and satellite waves, 
(e) (Main wave electrostatic energy/thermal energy) c 1,96 x 10 
(Test wave electrostatic energy/ thermal energy) ? 1,16 x 10 
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6.3 Theory 


At time t = 0 , the total electric field due to the main wave 
and the test wave is given by 

E(x,0) = {E o exp ik^x + exp i(k Q - A k ) x } + {c.c.} , (6.2) 

where E Q is the test wave amplitude, and c.c. denotes complex 
conjugate. When e is small, Eq. (6.2) can be written as 


E(x, 0) =* c 2 (x)E q exp[ik Q x - 6(x)] + {c.c.} , 

2 

c (x) = exp (e cos £kx) =- 1 + e cos A k x , (6.3) 

0(x) = e Sin £kx . 


This shows that the test wave can be regarded as (spatial) modulation 
of the amplitude and phase of the main wave when € « 1 . 

In the absence of modulation, the electric field of the main wave 
is given by 


S u (x,t) = E Q jexp “ J Q(t ' )dt exp[-i(uj Q t - k Q x)] + {c.c.}, 


(6.4) 


o(t') = - [i6u)(t ') + r(t')] , 


t' = V 1 


where fta) is the nonlinear frequency shift, and y is the damping rate 

of the main wave. It is assumed that gou and y are functions of 

amplitude and time only through the product mt , and that m is 

B B 

independent of time. The use of the form given by Eq . (6.4) would be 

Valid if 7 T /tW D « 1 , where y is the linear Landau damping rate of 
L B L 

48 

the main wave. 
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To incorporate the slow amplitude and phase variations in space 

2 

due to the presence of the test wave, is replaced by c (x)E^ , 

and the phase 6(x) defined in Eq. (6.3) is included. The modulated 
wave is then given by, 


c(x)<ju B t 

E n (x,t)_= c 2 (x)K 0 n<t')dt'] 

0 (6.5) 

X exp|-i[u> Q t - k Q X + 6(x)]j + {c.c.} . 

2 

Equation (6.5) gives a solution in space for a given amplitude, c (x)E Q , 
and initial phase, G(x) . If £k is given, the solution of Eq . (6.5) 
is correct only for time t < 2j\/^kv^ , where v^ is the phase 
velocity of the main wave. 

To obtain the temporal evolution of the Fourier modes, Eq. (6.5) is 
Fourier-transformed in space by the relation, 


;,t) =/c 


E (k, t ) = / dx exp(-ikx )E 

m ' 


We obtain, after some manipulation, 

|E in ( k ,t)| = E u (t) |6(k-k Q ) + e[ (l+A(t)J / 

r . 2 / . 2 . . -i 1 / 2 

+ € [A (t) + B (t)] ' 


m (x,t) . (6.6) 

+ B 2 (t)] 1 / 2 6(k-k o + &k) 
6(k-k Q -&k) + 0(e 2 ) J , 


(6.7) 


where 6( ) is the Dirac delta-f unction, and 
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Equation (6.7) shows that the main wave is unchanged to order € , 
and predicts the existence of a satellite at k Q + * Note that the 

satellite is linear in the test wave amplitude, eE^ , consistent with 
the results of the simulation shown in Fig. 6.1(a)- (c). The change in 
the growth rate to be seen in Fig. 6 . 1(c)- (e) when the main wave ampli- 
tude is varied, is suggested by Eq. (6.7), since A(t) and B(t) 
depend on the main wave evolution. 


6 . 4 Comparison with Simulation 

To make calculations from Eq . (6.7), we need to know values of 5oi 
and 7 to be substituted in Eq . ( 6 , 8 ). Although the nonlinear frequency 
shift, 6 u) , and growth rate, 7 , have been predicted theoretically , 48,66 
we prefer to use £u;(t) and 7 (t) determined from the results of Our 
simulation. This avoids error due to the observed main wave evolution 
not being exactly as these theories predict. To determine 6 u)(t) and 
7 (t) , we have first tabulated the phase change and amplitude of the 
complex Fourier amplitude of the main wave after every time- step. 
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Since these quantities contain large ripples, we have smoothed them 
using the least-square polynomial of degree one through five successive 
amplitudes. In Fig. 6.1 are plotted the theoretical calculations from 
Eq, (6.7) using these values. We see that there is good agreement 
between the theory and the simulation. In particular, the theoretical 
growth rate of the satellite at the earliest stage increases as the 
main wave amplitude increases, in very good agreement with the simulation. 
It is also significant that in each case the observed satellite energy 
level is in agreement with the calculated one. Although the calculated 
behavior of the test wave agrees well with the simulation, detailed 
observation shows that the simulated test wave first damps at the 
linear Landau damping rate, and then at the increased rate in agree- 
ment with the calculated one. 

6.5 Summary 

We have studied the temporal behavior of the satellite wave pro- 
duced when a large amplitude electron plasma wave and a small amplitude 
test wave are launched simultaneously. It has been shown that a theory 
which treats the test wave as a slow modulation of the amplitude and 
phase of the main wave explains well quantitatively the rapid growth 
and energy level of the satellite observed in our computer simulation. 
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7. CONCLUSIONS 


In Sections 2-6, we have studied linear and nonlinear phenomena 
associated with electron plasma waves, using a low- noise hybrid simu- 
lation model. Despite its attractive features, little use had previously 
been made of this model since it was proposed by Denavit in 1972. 

In Section 3, the model was first studied in detail, and demon- 
strated to simulate precisely the linear wave dispersion characteristics 
predicted by theory for long wavelength collective behavior. This 
verification of the validity and effectiveness of the simulation model 
is very important as a starting point for the subsequent study of non- 
linear phenomena. It also serves to establish the validity of the 
widely-used Cloud-in-Cell model, and the Langdon theory describing the 
finite-size particle model. Quantitative results in the very low 
energy range discussed here have never been obtained previously with 
such a high degree of accuracy with the simple particle models of 
Section 2. 

In Section 4, the low-noise model was used to investigate the 

behavior of a monochromat ic wave in both the linear and nonlinear 
* 

regimes. It was found that existing nonlinear theories are qualitatively 

in good agreement with the simulation results, but that there are some 

significant differences. In particular, the phase-mixing has been 

found to be slower than predicted. A new contribution of this section 

is a measurement of the nonlinear frequency shift, which is shown to 

1/2 

vary as E ' 

Section 5 was devoted to the investigation of the sideband, or 
trapped-particle, instability. Very good agreement was obtained for 
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CU B /Wp £ °* 2 between the results of the simulation, and a simple 

theory by Kruer et al. and a quasilinear theory by Brinca. We have not 
studied in detail the characteristics of the sideband instability for a 
heavily-damped main wave. This problem remains to be investigated 
further by both simulation and theory. 

In Section 6, we have studied nonlinear process involving coupling 
between a test wave and a large amplitude wave to produce a satellite 
wave. A simple theory, based on modulation of the large amplitude wave, 
was shown to explain the behavior of the satellite wave. This process 
may have important consequences in connection with the sideband insta- 
bility discussed in Section 5; when a test wave in one sideband is 
launched at a level above the fluctuation level, as is often done in 
experiments on the sideband instability, a corresponding wave in the 
other sideband may grow rapidly to a comparable level to that of the 
test wave, before the effects of the trapped particle instability 
discussed in Section 5 come into play. As a consequence, this nonlinear 
process may, for example, affect the measured energy spectrum indepen- 
dently of the sideband instability. 

We wish to emphasize in connection with Sections 3-6 that all of 
the simulations that have been presented were carried out under conditions 
for which the assumptions of relevant theoretical models could be 
approached, and in realistic energy ranges compared with those under 
which laboratory experiments are performed. It should be noted in 
relation to the latter, however, that our simulations have been 
concerned with temporal variations in a periodically bounded system, 
rather than with spatial variations in an effectively unbounded system. 
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Although we have investigated only one-dimensional problems, it 
seems straightforward to extend the hybrid simulation model to two and 
three dimensions. The effects of magnetic field could also be included 
at the cost of increased complication. The smoothing operation becomes 
more involved and time-consuming as the dimensionality is increased, 
and magnetic field is included. Even if it may not yet be economically 
feasible to extend its use to multidimensional problems with magnetic 
field included, the hybrid simulation model can serve very well, with 
reasonable cost, to achieve a very low fluctuation level given the 
capacity of currently available computers. 
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APPENDIX: DERIVATION OF WEIGHTING FUNCTIONS 


Derivation : 

0 

Consider a smoothing operation in velocity-space, 

f( v i ) f(v)w(v.-v) , (A . 1 ) 

V 

where denotes the i-th velocity grid point, and the summation is 

over particle velocity, v . The n-th order moment of the distribution 
function before smoothing is 


<v n > = y] v n f(v) . (A. 2) 

v 

After smoothing, the n-th order moment is given by 


<v n > v“ f(v.) 


(A. 3) 


Substituting Eq. (A.l) into Eq. (A. 3) and equating (v* 1 ) and (v* 1 ) 
yields 


w<v, - 


V) 


n 


v 


(A. 4) 


for any value of v . 

Velocity, v , can always be expressed as 

v = v . + fiv , 

J J 

= (j + r)£v (0 £ r s 1) , 

(A. 5) 
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where v . [= j&v] denotes the nearest velocity grid point such that 
J 

Vj £ v , JVj is the increment from the grid point, and is the 

velocity grid size. Substituting Eq . (A, 5) into Eq . (A. 4), and letting 
P = i - J , gives 

+ P) n w[(p - r)Av] = (j + r) n . (A. 6) 

P 

Using the binomial expansion, it will be seen that Eq. (A. 6) is satis- 
fied if 

E m m 

p w[(p - r)&v] = r (m = 0, 1, . . . n) . (A. 7) 

P 

Consider the Lagrangian interpolation of the function r™ , with 

89 

n + 1 points, given by 

s 

E m ( n+ 1 ) m / . _ . 

p A (r) = r , (A. 8) 

P=1~S 

where the A^ n+1 \r) are the Lagrangian coefficients with 0 £ r £ 1 , 
and s is an integer. Since m ^ n i the interpolation is exact. 
Comparing Eqs . (A. 7) and (A. 8), it follows that the desired weighting 
function may be given by 

w[(p - r) A v ] _ A^ n+1 \r) . (A. 9) 

Rewriting Eq . (A. 9) gives 

w(v) = A p n+1) (p - [<p - !)Av ^ V <: P Av] , (A. 10) 

where 1 - s ^ p ^ s . 
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When n is odd, the weighting function is given by Eq. (A. 10) 
with s = (n+l)/2 . When n is even, the Lagrangian coefficients do 
not give even weighting functions. In this case, they may be obtained 
by symmetrization as follows, 


w(v) = 


(n+i) 

/, n 

v \ 

1 r /n \ 

n 1 

-(n/2) 

( + 2 + 

Avj 

I [- \2 + 1 M v S v s - 2 Av j ’ 

A (n + 1)( 

>-b) 

+ A 

.‘"Wl-p + 

1-P \ A v/ 

j(p-l)Av ^ V s p^vj , 

<rn- 1 > / 

\ n 

v_) 

\ r n 1 

in \ *1 

A -(n/2)( 

, 1+ 2 - 

&V) 

1 [ 2 Av s v s ( 

,5 + *> Av J ■ 

(A. 11) 

e P £ S , 

and s 

= (n + 2)/2 . 



Although the weighting functions were derived for velocity-space, 
Eqs. (A. 10) and (A. 11) can be used for coordinate-space by replacing 
v by x . 


Examples : For n=l , Eq. (A. 10) is written as 

i a o 2) (" fe) = 1 + b (_Av s v 5 0) ' 


(D, ^ 

w (v) = 


'ff-b)- 


(A, 12) 


1 - 


&v 


(0 £ v £ £v) 


The smoothing operation using this weighting function conserves particles 
and momentum. 
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For nae2 , applying Eq. (A. 11) yields 


( 2 ), . 

w (v) 




1 



l - 


3v 

£\v 


- - ( — ) 
4 \ Av / 


(O ^ v £ Av) , 


— ) (Av <: v £ 2Av) . 
Av / 


(2 ) 

In the interval, -2£v £ v ^ 0 , w (v) is defined by symmetry* This 
weighting function conserves particles, momentum, and energy* The 
functions w^^(v) and w^ 2 ^(v) are shown in Pig. A,l. 



FIG. A . 1 . Linear (n = 1) and quadratic (n = 2) weighting 
functions. (Adapted from Fig. 3 of Ref. 6.). 
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